#include "mexUpperEnvelope.c"
void revert_ichoice(const int ichoice, int &de, int &dq, double &dk);

double interp_consumption(double optC_ss[NUM_coefEmax], double theta_new[], double logcprice_new, double qstock_new, int ia, int ieduc){
	
	double interp_c; 
		
	double state_new[DIM_smolyak]; 
	for (int isd=0; isd<DIM_smolyak; isd++){	
		state_new[isd] = get_stateSMOLYAK(isd, theta_new, logcprice_new, qstock_new, ia, ieduc); 
	}
	
	double tot_wc = 0.0; 
	double tot_w  = 0.0; 
	
	for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
		
		double state_ss[DIM_smolyak]; 								
		for (int idim=0; idim<DIM_smolyak; idim++){		
		state_ss[idim] = (GRID_smolyak[ipt][idim] + 1.0)*(MAX_coefEmaxstate[ia][ieduc][idim]-MIN_coefEmaxstate[ia][ieduc][idim])/2.0 + MIN_coefEmaxstate[ia][ieduc][idim];
		}
		
		
		double dis = 0.0; 
		for (int idim=0; idim<DIM_smolyak; idim++){
		
			dis += fabs (state_ss[idim] - state_new[idim]); 
		}
		
		
		if (dis < 1.0e-4){
		
			interp_c = optC_ss[ipt]; 
			
			break; 
		}
		
		
		double wt = 1.0/ sqrt( dis ); 
		
		#ifdef DEBUG
		if ( isnan(wt) || wt > 1.0e+5 || wt < 0.0 ) {
			printf("\n Errors in interp_consumption wt = %f, exiting \n", wt); exit(1);  
		} 
		#endif 
	
		tot_wc += wt * optC_ss[ipt]; 
		tot_w  += wt; 
		
		interp_c = tot_wc / tot_w; 
		
		
	} // ipt 
	
	
	#ifdef DEBUG
	if (interp_c <= 1000.0/1.2) {printf("\n interp_c %f < 1000 \n Error!!", interp_c); exit(1); }
	#endif
	
	return interp_c; 
}	


double get_appEGM_consumption(int de_prev0, int itype, int ieduccatpr0, int ia, int ieduc, double theta[], double logcprice, double qstock, double m, int ichoice){

	#ifdef DEBUG
	if (ia >= NUM_age || ia < 0 )      { printf("\n Error in get_appEGM_consumption: ia %2d \n ... exiting...", ia); exit(1); };		
	#endif
	
	int de_prev = de_prev0; int ieduccatpr = ieduccatpr0; 
	
		int n_de   = 2; 
		int educ_available = (ia < NUM_ageEduc ) * (ieduc < NUM_ieduc-1 ); 
		if (ia+INIT_age > MAX_ageHS && ieduc+MIN_educ < 12) educ_available = 0; 
		if (educ_available==0) n_de = 1; 				

		int n_ieduccatpr  = 1 * (ia >= NUM_age_ptr) + NUM_ieduccatpr * (ia < NUM_age_ptr); 
		
		if (n_de==1) {de_prev = 0;}
		if (n_ieduccatpr==1) {ieduccatpr=0; }
		
		
		double y_i[] = {0.0, 0.0}; int im_0=0; int im_1=0; // int ipt_0=0; int ipt_1 = 0; 
		double x_i [] = {0.0, 0.0}; 
		
		if (ia < NUM_ageEduc){
			
			int ia_ageYoung = ia; 
					
			find_GridIndex(m, Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr], NUM_MGrid, im_0, im_1); 
			
			x_i[0] = Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_0]; 
			x_i[1] = Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_1]; 
			
			y_i[0] = interp_consumption(Common_ageYoung_optC[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_0], theta, logcprice, qstock, ia, ieduc); 
			y_i[1] = interp_consumption(Common_ageYoung_optC[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_1], theta, logcprice, qstock, ia, ieduc); 
			
		} else {
			
			if (ichoice > 2){printf("\n Errors in get_appEGM_consumption \n exiting..\n"); exit(1);  }
			
			int ia_ageOld = ia - NUM_ageEduc; 
			
			find_GridIndex(m, Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype], NUM_MGrid, im_0, im_1); 
			
			x_i[0] = Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im_0]; 
			x_i[1] = Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im_1]; 
			
			y_i[0] = interp_consumption(  Common_ageOld_optC[ichoice][ia_ageOld][ieduc][itype][im_0], theta, logcprice, qstock, ia, ieduc); 
			y_i[1] = interp_consumption(  Common_ageOld_optC[ichoice][ia_ageOld][ieduc][itype][im_1], theta, logcprice, qstock, ia, ieduc); 
			
		}
					
		double y = y_i[0] + ( m - x_i[0] ) * (y_i[1] - y_i[0])/(x_i[1] - x_i[0]); 
		
		
		int de, dq; double dk; 
		revert_ichoice(ichoice, de, dq, dk); 
		int educ = ieduc + MIN_educ; 		
		double probs = get_survivalProb(ia, educ, qstock, dq, de); 
		double delta = get_delta(theta, qstock); 
		
		y *= pow( probs * delta, -1.0/cu_g ); 
		
		
	#ifdef DEBUG

		double state[DIM_smolyak]; 
		
		for (int isd=0; isd<DIM_smolyak; isd++){
		
			state[isd] = get_stateSMOLYAK(isd, theta, logcprice, qstock, ia, ieduc); 
		}
						
		for (int idim=0; idim<DIM_smolyak; idim++){
			
			if ( isnan(state[idim]) || state[idim] < MIN_coefEmaxstate[ia][ieduc][idim] - 1.0e-4 || state[idim] > MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-4 ){		
				printf("\n\n Errors in get_appEGM_consumption(%2d): idim %2d, state %f, [%f, %f] \n\t ieduc %2d, theta [%f, %f], m %f, maxdebt %.2f \n ... exiting ... \n", ia,
				idim, state[idim], MIN_coefEmaxstate[ia][ieduc][idim], MAX_coefEmaxstate[ia][ieduc][idim], ieduc, theta[0], theta[1], m, MAX_debt[ia][ieduc] ); 
				exit(1); 
			}
			
			if (state[idim] > MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-10) state[idim] = MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-10; 

		}
		

		if ( isnan(y) || y_i[0] < 0.0 || y_i[1] < 0.0 || ia > NUM_age-1){
			
			printf("\n\n Errors in get_appEGM_consumption: ichoice %d, age %2d, ia %2d , ieduc %d, itype %d, y_i [%f, %f], x_i [%f, %f], state [%f %f %f %f], im [%d, %d]" , ichoice, ia+INIT_age, ia, ieduc, itype, y_i[0], y_i[1], x_i[0], x_i[1], state[0], state[1], state[2], state[3], im_0, im_1 ); 

		
			exit(1); 	

		}
	
	#endif 
		
	return y; 

		
}		



double get_appEGM_v(double coef_ageYoung[NumChoice_ageYoung][NUM_ageEduc][NUM_de][NUM_ieduc][NUM_isavingspr][NUM_ieduccatpr][NUM_MGrid][NUM_coefEmax], double coef_ageOld[NumChoice_ageOld][NUM_ageNoEduc][NUM_ieduc][NUM_isavingspr][NUM_MGrid][NUM_coefEmax], int de_prev0, int itype, int ieduccatpr0, int ia, int ieduc, double theta[], double logcprice, double qstock, double m, int ichoice){

	#ifdef DEBUG
	if (ia >= NUM_age )      { printf("\n Error in get_appEGM_v(): ia %2d \n ... exiting...", ia); exit(1); };		
	#endif
	
	int de_prev = de_prev0; int ieduccatpr = ieduccatpr0; 
	
		int n_de   = 2; 
		int educ_available = (ia < NUM_ageEduc ) * (ieduc < NUM_ieduc-1 ); 
		if (ia+INIT_age > MAX_ageHS && ieduc+MIN_educ < 12) educ_available = 0; 
		if (educ_available==0) n_de = 1; 				

		int n_ieduccatpr  = 1 * (ia >= NUM_age_ptr) + NUM_ieduccatpr * (ia < NUM_age_ptr); 
		
		if (n_de==1) {de_prev = 0;}
		if (n_ieduccatpr==1) {ieduccatpr=0; }
		

		double state[DIM_smolyak]; 
		
		for (int isd=0; isd<DIM_smolyak; isd++){
		
			state[isd] = get_stateSMOLYAK(isd, theta, logcprice, qstock, ia, ieduc); 
		}
						
		for (int idim=0; idim<DIM_smolyak; idim++){
			
			if ( isnan(state[idim]) || state[idim] < MIN_coefEmaxstate[ia][ieduc][idim] - 1.0e-4 || state[idim] > MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-4 ){		
				printf("\n\n Errors in get_appEGM_v(%2d): idim %2d, state %f, [%f, %f] \n\t ieduc %2d, theta [%f, %f], m %f, maxdebt %.2f \n ... exiting ... \n", ia,
				idim, state[idim], MIN_coefEmaxstate[ia][ieduc][idim], MAX_coefEmaxstate[ia][ieduc][idim], ieduc, theta[0], theta[1], m, MAX_debt[ia][ieduc] ); 
				exit(1); 
			}
		

			if (state[idim] > MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-10) state[idim] = MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-10; 
			
		}
		
		struct_smolyak smolyak; 
		
		double y_i[] = {0.0, 0.0}; int im_0=0; int im_1=0; // int ipt_0=0; int ipt_1 = 0; 
		double x_i [] = {0.0, 0.0}; 
		
		 
		if (ia < NUM_ageEduc){
			
			int ia_ageYoung = ia; 
					
			find_GridIndex(m, Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr], NUM_MGrid, im_0, im_1); 
			
			x_i[0] = Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_0]; 
			x_i[1] = Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_1]; 
			
			y_i[0] = smolyak.PolyInterp(ChebyIndex, state, MIN_coefEmaxstate[ia][ieduc], MAX_coefEmaxstate[ia][ieduc], coef_ageYoung[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_0]);
			y_i[1] = smolyak.PolyInterp(ChebyIndex, state, MIN_coefEmaxstate[ia][ieduc], MAX_coefEmaxstate[ia][ieduc], coef_ageYoung[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im_1]);
			
			
		} else {
			
			if (ichoice > 2){printf("\n Errors in get_appEGM_v \n exiting..\n"); exit(1);  }
			
			int ia_ageOld = ia - NUM_ageEduc; 
			
			find_GridIndex(m, Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype], NUM_MGrid, im_0, im_1); 
			
			x_i[0] = Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im_0]; 
			x_i[1] = Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im_1]; 
			
			y_i[0] = smolyak.PolyInterp(ChebyIndex, state, MIN_coefEmaxstate[ia][ieduc], MAX_coefEmaxstate[ia][ieduc], coef_ageOld[ichoice][ia_ageOld][ieduc][itype][im_0]);
			y_i[1] = smolyak.PolyInterp(ChebyIndex, state, MIN_coefEmaxstate[ia][ieduc], MAX_coefEmaxstate[ia][ieduc], coef_ageOld[ichoice][ia_ageOld][ieduc][itype][im_0]);
			
			
		}
			
		double y = y_i[0] + ( m - x_i[0] ) * (y_i[1] - y_i[0])/(x_i[1] - x_i[0]); 
		
	#ifdef DEBUG

		if ( isnan(y) || ia > NUM_age-1){

			printf("\n\n Errors in get_appEGM_v(): ichoice %d, age %2d, ia %2d , ieduc %d, itype %d, y_i [%f, %f], x_i [%f, %f], state [%f %f %f %f], im [%d, %d]" , ichoice, ia+INIT_age, ia, ieduc, itype, y_i[0], y_i[1], x_i[0], x_i[1], state[0], state[1], state[2], state[3], im_0, im_1 ); 
			
		
			exit(1); 	

		}
	
	#endif 
		
	return y; 

		
}		
		

void revert_ichoice(const int ichoice, int &de, int &dq, double &dk){
	
	switch(ichoice) {
	
	case 0: 
		de = 0; 
		dq = 0; 
		dk = 1.0; 
		break; 
		
	case 1: 
		de = 0; 
		dq = 1; 
		dk = 1.0; 
		break; 
	
	case 2: 
		de = 1; 
		dq = 0; 
		dk = 0; 
		break; 
		
	case 3: 
		de = 1; 
		dq = 0; 
		dk = 0.5; 
		break; 
	
	case 4: 
		de = 1; 
		dq = 1; 
		dk = 0.0; 
		break; 
	
	case 5: 
		de = 1; 
		dq = 1; 
		dk = 0.5; 
		break; 
	
	default : 
		printf("\n Error in revert_ichoice(), ichoice %d ", ichoice); exit(1); 					
		break; 
	
	}

}
int  get_ichoice(const int de, const int dq, const double dk){
	
	int ichoice = -9;
	 
	if (de==0 && dq ==0){
	ichoice = 0; 
	} else if (de==0 && dq ==1){
	ichoice = 1; 
	} else if (de==1 && dq ==0 && dk < 0.4){
	ichoice = 2; 
	} else if (de==1 && dq ==0 && dk > 0.4 && dk < 0.9){
	ichoice = 3;  
	} else if (de==1 && dq ==1 && dk < 0.4){
	ichoice = 4; 
	} else if (de==1 && dq ==1 && dk > 0.4 && dk < 0.9){
	ichoice = 5; 
	} 
	
	#ifdef DEBUG
	if (ichoice==-9){ printf("\n Errors in get_ichoice: de %d, dq %d, dk %f \n ", de, dq, dk); exit(1); }
	#endif 
	
	return ichoice; 

}



inline int getGovSubRule(int i_model, int de, int dq, double theta[], int educ, int ia){
	
	int isub = 0; 	
	
	if (cc_taxPolicy == 1) {
	
	if (i_model == 0){
		
		isub = (de==1); 
		
	} else if ( i_model == 1 ){
		
		isub = (theta[0] < 0.0)*(theta[1] < 0.0)*(ia < NUM_ageEduc); 
	
	} else if ( i_model == 2 ){
		
		isub = (dq==0) * (ia < NUM_ageEduc);
		  
	} else {
		
		printf("\n\n Error in getGovSubRule: i_model = %d \n exiting...\n", i_model);
	}
	
	}
	
	return isub; 
	
}


double get_totalwealthEGM(const int de, const double dk, const int dq, const struct_infoset &infoset, double &gtrcmin){
	
	int    educ = infoset.educ; 
	double s    = infoset.savings; 
    double hrs = 0.0 + HOURSP * (dk>0.1 && dk < 0.9) + HOURS * (dk > 0.9);
	double inc_wages  = get_wagerate(dk, de, infoset) * hrs; 
	

	double inc_s = R_lend * s * (s>0.0) + R_borrow * s * (s<0.0) ; 

	if (educ >= 12 && de ==1 && s < 0.0) { 
	    int ieduc = educ-MIN_educ;
		double debt_gsl = min(-s, MAX_debt_gsl[ieduc] ); 
		
		if (educ < 16 ){
			inc_s += R_borrow * debt_gsl; 
		}
		    
	}
	
	
	double c_pt = getc_subsidy(infoset.ia, educ, de); 
	double ptr = get_ptransfer(de, dk, infoset);
	
	
	double y = (c_pt + ptr) + inc_wages + ( s + inc_s) ;  	
	


	gtrcmin = 0.0; 
	if (de==0 && s <= 0.0 && (c_pt + ptr) + inc_wages < cu_min) gtrcmin = cu_min - (c_pt + ptr + inc_wages); 
	

	
	double cost_educ = 0.0; 
	if (de==1 && educ >=12) {
		int itc = educ+de - 13; 
		cost_educ = (cbc_tc[itc] - cbc_grants);
	}
	
	double cost_dq = infoset.health; 
	
	double subtax = 0.0; 
	if (cc_taxPolicy ==1){
			
			subtax =  cbc_sub * getGovSubRule(0, de, dq, infoset.theta, educ, infoset.ia) 
				+ gtr_sub * getGovSubRule(1, de, dq, infoset.theta, educ, infoset.ia) 
				+ gtr_csub * getGovSubRule(2, de, dq, infoset.theta, educ, infoset.ia) 
				- dq * cost_dq * TAXrate_dq - get_incomeTax(inc_wages, inc_s); 
	
	}
	
	// get m from savings and choices 
	return ( gtrcmin + y - cost_educ*de - cost_dq * dq + subtax ); 
	
}

double get_Croom(int ia, int educ, int de){
    double roomboard = 0.0;
    if (de == 1 && educ+de>=13) {
        int itc = educ+de - 13; roomboard = cbc_roomboard[itc];
    }
    return roomboard;  
}


double get_minS_EGM(int ia, int educ, int de, double dk, const double debtlimit_p){
//--- This function returns a lower bound of savings depends on the individual's state
	
	double debt_g = 0.0; 
	
	int itc = educ + de - 13; 
	if (de==1 && itc >=0) {
        
        for (int i =0; i<= itc; i++){
        	debt_g += min( max_debtg_e[itc], cbc_tc[itc] - cbc_grants + cbc_roomboard[itc] );
        }
		
		double max_debtg1 = max_debtg[0]+(educ+de>16)*max_debtg[1]; 
		if (debt_g > max_debtg1 ) debt_g = max_debtg1; 
		
	}
	
	if (de==0 && educ>=16 && ia + INIT_age <= 25 ) {
		
		debt_g += maxpost_debtg; 
		
	} 
	
	#ifdef DEBUG
	if (itc >=8 ) { printf("\n\n Errors in tuition cost and debt ... exiting ... \n\n"); exit(1); }
	#endif 
		    
	double min_s = min( -1.0 * debt_g, -1.0 * debtlimit_p); 
		
	#ifdef DEBUG 
		if ( debtlimit_p < -1.0e-4 || min_s > 1.0e-4 ) { printf("\n\n Errors in get_minS: %f , debt_p %f \n \n exiting \n\n", min_s, debtlimit_p ); exit(1); } 
	#endif 
	
	return ( min_s ); 
	
}



void Interp_ValuePolicy_EGM_new(double &m, double &v, double &c, double &snext, double &dv, double &gtransfer, double &gtrcmin, double &debtlimit_p, double &sp_lo, double &consumption_min, const int de, const double dk, const int dq, const struct_infoset &infoset){

	// initialization 
	m = -1.0e+30; 
	v = -1.0e+30; 
	c = -1.0e+30; 
	dv = 0.0; 
	gtransfer = 0.0; gtrcmin = 0.0; debtlimit_p = 0.0; 


	int ia = infoset.ia; 
	int ieduc = infoset.educ-MIN_educ; 
//	const int 	 ieduc_next  = ieduc + de; 
//	const int 	 educ_next   = ieduc_next + MIN_educ; 
//	const double kstock_next = get_knext(infoset.ia, infoset.educ, infoset.kstock, dk, de);			
	
	
	double state[DIM_smolyak]; 
	for (int isd=0; isd<DIM_smolyak; isd++){
		double logh = infoset.logcprice; 
		#ifdef INFOSET_hnext
		logh = get_logcprice_next(infoset.logcprice, infoset.eps_health_next); // logcprice_next
		#endif 
		state[isd] = get_stateSMOLYAK(isd, infoset.theta, logh, infoset.qstock, ia, ieduc); 
	}	
	for (int idim=0; idim<DIM_smolyak; idim++){
		if ( isnan(state[idim]) || state[idim] < MIN_coefEmaxstate[ia][ieduc][idim] - 1.0e-4 || state[idim] > MAX_coefEmaxstate[ia][ieduc][idim] + 1.0e-4 ){		
		printf("\n\n Errors in get_appEmax(%2d): idim %2d, state %f, [%f, %f] \n\t ieduc %2d, theta [%f, %f], savings %f, maxdebt %.2f \n ... exiting ... \n", ia,
		idim, state[idim], MIN_coefEmaxstate[ia][ieduc][idim], MAX_coefEmaxstate[ia][ieduc][idim], ieduc, infoset.theta[0], infoset.theta[1], infoset.savings, MAX_debt[ia][ieduc] ); 
		exit(1); 
		}
	}
	
	
    debtlimit_p  = get_debtlimit_ptotal(infoset.theta, ia, infoset.educ, de, dq, infoset.qstock); 
	sp_lo = get_minS_EGM(ia, infoset.educ, de, dk, debtlimit_p); 
	m = get_totalwealthEGM(de, dk, dq, infoset, gtrcmin) - sp_lo;      // includes consumption subsidy 	
	
	double c_ptsubsidy = getc_subsidy(ia, infoset.educ, de); 
	double c_minexpenditure = getc_minexpenditure(ia,  infoset.educ, de, c_ptsubsidy);
	consumption_min  = max(c_ptsubsidy, c_minexpenditure); 
	
	
	int ichoice = get_ichoice(de, dq, dk); 
	
	double interest_rate = (1.0 + R_lend*(infoset.savings>=0.0) + R_borrow*(infoset.savings<0.0)); 
		
	// choice is available if there is sufficient resources 
	if ( m >= consumption_min ){

	
		c  = get_appEGM_consumption(infoset.de_prev, infoset.isavingspr, infoset.ieduccatpr, infoset.ia, ieduc, infoset.theta, infoset.logcprice, infoset.qstock, m, ichoice); 
		
		if ( c < consumption_min ) c = consumption_min; 
		
		snext = m + sp_lo - c; 
		
		if (snext < sp_lo) { snext = sp_lo;  c = m + sp_lo - snext; }
		
		v  = - dq * infoset.eps_uq * sigeps_uq - de * infoset.eps_ue * sigeps_ue + get_appEGM_v(CoefV_ageYoung_EGM, CoefV_ageOld_EGM, infoset.de_prev, infoset.isavingspr, infoset.ieduccatpr, infoset.ia, ieduc, infoset.theta, infoset.logcprice, infoset.qstock, m, ichoice); 	
		
		dv = interest_rate * get_dUdC_EGM(c, ia, infoset.educ); 
	
	} else if ( m - c_ptsubsidy < 0.0 && dq == 1 ) { // not enough money to buy
      
        m = -1.0e+30;
        v = -1.0e+30;
        c = -1.0e+30;
        dv = 0.0;
        gtransfer = 0.0; gtrcmin = 0.0;
        
    } else if (de==0  || (de==1 && c_ptsubsidy >= c_minexpenditure) ){ // can happen if youth is heavily in debt 
				
		c  = max(m, c_minexpenditure); 
		
		snext = m + sp_lo - c; 
		
		// readjust c according to consumption subsidy		
		if ( c < c_ptsubsidy ) c = c_ptsubsidy; 
		
		v  = - dq * infoset.eps_uq * sigeps_uq - de * infoset.eps_ue * sigeps_ue + get_appEGM_v(CoefV_ageYoung_EGM, CoefV_ageOld_EGM, infoset.de_prev, infoset.isavingspr, infoset.ieduccatpr, infoset.ia, ieduc, infoset.theta, infoset.logcprice, infoset.qstock, m, ichoice); 	
				
		dv = interest_rate * get_dUdC_EGM(c, ia, infoset.educ); 
		 		
	} 
	
    gtransfer = gtrcmin; 
	
}	


double get_EvnextEGM_new(double theta[], int ia, int ieduccatpr0, int isavingspr, int de_prev0, int ieduc, double kstock0, double savings, double logcprice_foresee, double qstock, double &emax_dv) {	
// Note: logcprice_foresee = E[ log health(t) | t-1 ]

	// output variable initialization 
	double emaxv = 0.0;  
	emax_dv = 0.0; 
	
	// !!!!!!!!!!! Start: Defining deterministic portion of the information set !!!!!!!!!!!!!!
	struct_infoset infoset;
		
	int educ = ieduc + MIN_educ; 
	int ieduc5 = 1*(educ==12) + 2 * (educ >=13) * (educ<=15) + 3 * (educ==16) + 4 * (educ>16); 
			
	int de_prev = de_prev0; int ieduccatpr = ieduccatpr0; 
	
	int n_de   = 2; 
	int educ_available = (ia < NUM_ageEduc ) * (ieduc < NUM_ieduc-1 ); 
	if (ia+INIT_age > MAX_ageHS && ieduc+MIN_educ < 12) educ_available = 0; 
	if (educ_available==0) n_de = 1; 				

	int n_ieduccatpr  = 1 * (ia >= NUM_age_ptr) + NUM_ieduccatpr * (ia < NUM_age_ptr); 
	
	if (n_de==1) {de_prev = 0;}
	if (n_ieduccatpr==1) {ieduccatpr=0; }
		
	double kstock = get_kstock(educ, ia); 				
	double dqprev = (qstock>=1); double qstockprev = 0.0;  // approximate 
	double debtlimitLag = get_debtlimit_ptotal(theta, ia-1, educ-de_prev0, de_prev0, dqprev, qstockprev); 
	
										 	   				    	
	//--- Deterministic components except health, logcprice, delta, and skill 
	infoset.theta   = theta; 					 													
	infoset.ia = ia; 
	infoset.de_prev = de_prev;
	infoset.debtlimitLag = debtlimitLag;			

	infoset.educ    = educ; 
	infoset.kstock  = kstock;
	infoset.qstock  = qstock;  
											
	infoset.savings = savings; 
					
	infoset.delta   = get_delta(infoset.theta, infoset.qstock); 						

	infoset.ieduccatpr = ieduccatpr; 
	
	infoset.isavingspr = isavingspr; 
	
	// !!!!!!!!!!! Complete: Defining deterministic portion of the information set !!!!!!!!!!!!!!	
	int n_eps_health = neps_health; 
	
	double w_eps_health[n_eps_health], x_eps_health[n_eps_health]; 
	if (n_eps_health==1){
			w_eps_health[0] = 1.0; 
			x_eps_health[0] = 0.0; 
			
	} else {
		for (int ih =0; ih < neps_health; ih ++) {
			w_eps_health[ih] = weps_health[ih]; 
			x_eps_health[ih] = xeps_health[ih]; 
		}
	
	}
	
	
	for (int ih=0; ih < n_eps_health; ih++){
	for (int ie=0; ie < NUM_EPS_earnings; ie++){ 

	
		double wgt = WEPS_logwf[ieduc5][ie]  *  w_eps_health[ih];
		
		infoset.eps_logwf        = XEPS_logwf[ieduc5][ie];  
		
		
		//--- current price 
		double logcprice = logcprice_foresee + x_eps_health[ih] * sigeps_h;
		 
		if ( logcprice < MIN_logcprice ) {
		logcprice = MIN_logcprice; 
		} else if (logcprice > MAX_logcprice){
		logcprice = MAX_logcprice; 
		}
		infoset.eps_health_next = 0.0; 										
		
		
		// next period's price is known this period via eps_health_next
		#ifdef INFOSET_hnext	 
				
		logcprice = logcprice_foresee;
		
		double logcprice_next = get_logcprice_next(logcprice, x_eps_health[ih]); 
		
		infoset.eps_health_next  = x_eps_health[ih];
		
		#endif 
		
								
		infoset.health  = exp(logcprice); 	
		infoset.logcprice  = logcprice; 
		
		
		//--- skill 			
		infoset.skill           = get_skill(educ, theta, ia, infoset.eps_logwf); 
		
		
		int issp_eh = ie*ih; 
		infoset.eps_ptr   		= simp_rng_eps[issp_eh][8];  
		infoset.eps_ptrdummy   	= simp_rng_eps[issp_eh][9];  		
		#ifdef DEBUG
		if ( issp_eh > NUM_simsp ) { printf("\n\n Errors in get_simEmax(): issp_eh %d > NUM_simsp %d \n exiting...\n", issp_eh, NUM_simsp); exit(1); } 
		#endif 
		
		
		// assign preference shock to be zero
		infoset.eps_ue 	 		= 0.0;  					
		infoset.eps_uq   		= 0.0;  
		
				
		int educ_available = (ia < NUM_ageEduc ) * (ieduc < NUM_ieduc-1 ); 
		if (ia+INIT_age > MAX_ageHS && ieduc+MIN_educ < 12) educ_available = 0; 
		
		int nchoices = 2 + 4 * (educ_available==1); 
		
		// the output given the realization of shocks
		double maxV = 0.0; double dmaxV = 0.0; 
		
		// available choices based on model restriction. 
		double vvec[nchoices], duvec[nchoices]; // mvec[nchoices], cvec[nchoices], 
		for (int ichoice  = 0; ichoice  < nchoices; ichoice ++){			
			
			int de, dq; double dk; 
			revert_ichoice(ichoice, de, dq, dk); 	
			
			double m, v, c, snext, dv, gtransfer, gtrcmin, debtlimit_p, sp_lo, consumption_min; 
			Interp_ValuePolicy_EGM_new(m, v, c, snext, dv, gtransfer, gtrcmin, debtlimit_p, sp_lo, consumption_min, de, dk, dq, infoset); 
			
			vvec[ichoice]    = v; 
		//	cvec[ichoice]    = c; 
			duvec[ichoice]   = dv; 
		}
		
		
			double v_de0_dq0 = vvec[0]; double v_de0_dq1 = vvec[1];
			double dv_de0_dq0 = duvec[0]; double dv_de0_dq1 = duvec[1];
			
			
			double maxVe0  = v_de0_dq0; 
			double dmaxVe0 = dv_de0_dq0; 
			
		//	if (maxVe0<0.0) { maxVe0 = 0.0; }
			
		//	if (v_de0_dq1>0.0)
			{
			// ... 1. get maxVe0 = max (  v_de0_dq0, v_de0_dq1 - sigeps_uq) 
			//                   = v_de0_dq0 + ( v_de0_dq1-v_de0_dq0 - sigeps_uq ) * 1( cutoff - sigeps_uq > 0), where cutoff = v_de0_dq1-v_de0_dq0. 
			maxVe0 = get_maxV_2choices(sigeps_uq, v_de0_dq0, v_de0_dq1, dv_de0_dq0, dv_de0_dq1, dmaxVe0); 
			}			
						
			// ... 2. initialization 
			maxV  = maxVe0 ; 
			dmaxV = dmaxVe0 ; 
			
			
			
			if (nchoices == 6){
				
						double maxVe1, dmaxVe1;
				
						double v_de1_dq0 = -1.0e+30; double dv_de1_dq0 = -999.0; 
						double v_de1_dq1 = -1.0e+30; double dv_de1_dq1 = -999.0; 

						
						// de=1, dq=0
						if (vvec[2] >= vvec[3]){
							v_de1_dq0  = vvec[2];
							dv_de1_dq0 = duvec[2];  
						} else {
							v_de1_dq0  = vvec[3];
							dv_de1_dq0 = duvec[3];  
						}
						
						// de=1, dq=1
						if (vvec[4] >= vvec[5]){
							v_de1_dq1  = vvec[4]; 
							dv_de1_dq1 = duvec[4];
						} else {
							v_de1_dq1  = vvec[5]; 
							dv_de1_dq1 = duvec[5];	
						}
						
						maxVe1 = v_de1_dq0; 
						
					//	if (v_de1_dq1>0.0)
						{
						// ... 3. get maxVe1 = max (  v_de1_dq0, v_de1_dq1 - sigeps_uq); 
						maxVe1 = get_maxV_2choices(sigeps_uq, v_de1_dq0, v_de1_dq1, dv_de1_dq0, dv_de1_dq1, dmaxVe1); 
						}
						
					//	if (maxVe1>0.0)
						{
						// ... 4. update maxV = max (  maxVe0, maxVe1 - sigeps_ue); 
						maxV = get_maxV_2choices(sigeps_ue, maxVe0, maxVe1, dmaxVe0, dmaxVe1, dmaxV); 
						
						}
							
			} // educ_available == 1	
		
				
		emaxv  += ( wgt * maxV );  
		emax_dv += ( wgt * dmaxV ); 
				
	} // cigarrate price shock 
	}// earnings shock
	
	
	// || emaxv < 1.0e-30
	
	#ifdef DEBUG
	if (emax_dv <  1.0e-30 ){
		printf("\n Errors in get_EvnextEGM(), emaxv %f, emaxdu %f ", emaxv, emax_dv); 
	}
	#endif 
	
	return emaxv; 
		
}


double get_Sp_foc(double snext, void *params){
// 		This is the first order condition: lambda = dU/dC - delta* E(dV+/dS+)
// 		Notice: Sp exist if and only if dEmax_dSp > 0 
		
		struct struct_foc *p = (struct struct_foc *) params;
		
		double cashathand = p-> cashathand; 

		int de 	  = p-> de; 
		double dk = p-> dk; 
		int    dq = p-> dq; 
		
		struct_infoset infoset = p->infoset; 
		
		const int ia      = infoset.ia; 
		const int itype   = infoset.isavingspr; 
		double logcprice  = infoset.logcprice; 
		int    educ_next   = infoset.educ + de; 

		double c     = cashathand - snext;


		const int 	 ieduc_next  = educ_next - MIN_educ; // ieduc + de; 
		const double kstock_next = get_knext(infoset.ia, infoset.educ, infoset.kstock, dk, de);			
		
		double qstock_next = get_qnext_itype(infoset.qstock, dq, infoset.ia, itype);

		double loghnext_expected = logcprice; 
		{ 	
			loghnext_expected = get_logcprice_next(logcprice, 0.0); // infoset.eps_h
		} 
		
    	double probs = get_survivalProb(ia, infoset.educ, infoset.qstock, dq, de); 
		double delta = get_delta(infoset.theta, infoset.qstock); // infoset.delta; 
		
		double Edv_next  = -9.0;  	double ev_next = -9.0; 
	
		if (ia < NUM_age-1){ 
		
		ev_next = get_EvnextEGM_new(infoset.theta, ia+1, infoset.ieduccatpr, itype, de, ieduc_next, kstock_next, snext, loghnext_expected, qstock_next, Edv_next); 	 
		
		} else {
	
		ev_next  = get_PDV_EGM_new(ia+1, ieduc_next, infoset.theta, kstock_next, exp(loghnext_expected), snext, itype, dq, qstock_next, Edv_next); 
	
		}
		
		// Euler equation should satisfy: 
		
		double Inv_dVnext = inv_dUdC_EGM( Edv_next, infoset.ia, infoset.educ);  double c_new = pow( probs * delta, -1.0/cu_g ) * Inv_dVnext; 

		double lambda = c_new - c; 
		
		return lambda; 
}


double get_valueEGM(double snext, double c, int de, double dk, int dq, const struct_infoset &infoset, double &util, double &ev_next, double &probs, double &delta){
	
	const int ia      = infoset.ia; 
	const int educ    = infoset.educ; 	
	const int ieduc   = educ - MIN_educ;
	const int itype   = infoset.isavingspr; 
	const int ieduccatpr = infoset.ieduccatpr; 
	double logcprice  = infoset.logcprice; 
	
		const int 	 ieduc_next  = ieduc + de; 
		const double kstock_next = get_knext(infoset.ia, infoset.educ, infoset.kstock, dk, de);			
		
		double qstock_next = get_qnext_itype(infoset.qstock, dq, infoset.ia, itype);

		double loghnext_expected = logcprice; 
		{ 	
			loghnext_expected = get_logcprice_next(logcprice, 0.0); // infoset.eps_h
		} 
		
    	probs = get_survivalProb(ia, educ, infoset.qstock, dq, de); 
		delta = get_delta(infoset.theta, infoset.qstock); // infoset.delta; 

		double Edv_next  = -9.0;  	// double ev_next = -9.0; 
	
		if (ia < NUM_age-1){ 
		
		ev_next = get_EvnextEGM_new(infoset.theta, ia+1, ieduccatpr, itype, de, ieduc_next, kstock_next, snext, loghnext_expected, qstock_next, Edv_next); 	 
		
		} else {
	
		ev_next  = get_PDV_EGM_new(ia+1, ieduc_next, infoset.theta, kstock_next, exp(loghnext_expected), snext, itype, dq, qstock_next, Edv_next); 
	
		}
		
		// Euler equation should satisfy: double Inv_dV = inv_dUdC_EGM( EdV_next, ia, educ);  double c = pow( probs * delta, -1.0/cu_g ) * Inv_dV; 
		util = get_utilityEGM(de, dk, c, dq, ia, infoset.de_prev, ieduc, ieduccatpr, infoset.theta, infoset.qstock, infoset.eps_ue, infoset.eps_uq); 
		double v = util + probs * delta * ev_next; 

		return v; 

}

double get_valueEGM(double snext, double c, int de, double dk, int dq, const struct_infoset &infoset){
	
	const int ia      = infoset.ia; 
	const int educ    = infoset.educ; 	
	const int ieduc   = educ - MIN_educ;
	const int itype   = infoset.isavingspr; 
	const int ieduccatpr = infoset.ieduccatpr; 
	double logcprice  = infoset.logcprice; 
	
		const int 	 ieduc_next  = ieduc + de; 
		const double kstock_next = get_knext(infoset.ia, infoset.educ, infoset.kstock, dk, de);			
		
		double qstock_next = get_qnext_itype(infoset.qstock, dq, infoset.ia, itype);

		double loghnext_expected = logcprice; 
		{ 	
			loghnext_expected = get_logcprice_next(logcprice, 0.0); // infoset.eps_h
		} 
		
    	double probs = get_survivalProb(ia, educ, infoset.qstock, dq, de); 
		double delta = get_delta(infoset.theta, infoset.qstock); // infoset.delta; 

		double Edv_next  = -9.0;  	double ev_next = -9.0; 
	
		if (ia < NUM_age-1){ 
		
		ev_next = get_EvnextEGM_new(infoset.theta, ia+1, ieduccatpr, itype, de, ieduc_next, kstock_next, snext, loghnext_expected, qstock_next, Edv_next); 	 
		
		} else {
	
		ev_next  = get_PDV_EGM_new(ia+1, ieduc_next, infoset.theta, kstock_next, exp(loghnext_expected), snext, itype, dq, qstock_next, Edv_next); 
	
		}
		
		// Euler equation should satisfy: double Inv_dV = inv_dUdC_EGM( EdV_next, ia, educ);  double c = pow( probs * delta, -1.0/cu_g ) * Inv_dV; 
		
		double v = get_utilityEGM(de, dk, c, dq, ia, infoset.de_prev, ieduc, ieduccatpr, infoset.theta, infoset.qstock, infoset.eps_ue, infoset.eps_uq) + probs * delta * ev_next; 

		return v; 

}

/*
double solve_optSp(double &v_best, double &foc_xr, const double x0_lo, const double x0_hi, double cashathand, int de, double dk, const int dq, const struct_infoset &infoset){
// Solve for optimal Cn and Sp for given Mp, Np, K and rest_inc. 
	
	struct_foc params; 
	
	params.cashathand = cashathand; 

	params.isavingspr_next   = infoset.isavingspr;
	params.ieduccatpr_next   = infoset.ieduccatpr;
	
	params.infoset = infoset; 

	params.de   	= de; 	
	params.dq   	= dq; 	
	params.dk   	= dk; 	
	
	
	
	double foc_lo, foc_hi, v0_lo, v0_hi; int isbad; 
	
	
	
    double x_r = x0_hi;
		
    int status, status0,iter, max_iter;
    
//--foc (lambda) is increasing in snext: foc_lo < 0, foc_hi > 0.0
	double x_lo = x0_lo; 
	double x_hi = x0_hi; 	
	foc_lo = get_Sp_foc(x_lo, &params); 
	foc_hi = get_Sp_foc(x_hi, &params);  




    v0_lo   = get_valueEGM(x0_lo, params.cashathand-x0_lo, params.de, params.dk, params.dq, params.infoset);    
    v0_hi   = get_valueEGM(x0_hi, params.cashathand-x0_hi, params.de, params.dk, params.dq, params.infoset);
    
    // ... isbad = 0 (good interior solution), isbad = 1 (lower bound corner solution), isbad = 2 (higher bound corner solution), isbad = 3 (foc is wrong; undertermined)
    isbad = 3;
    
	#ifdef DEBUG_Sp_foc
          	printf("\n\t range [%f, %f], foc [%f, %f]\n", x_lo, x_hi, foc_lo, foc_hi);
    #endif
           		
	if (foc_hi > 0.0 && foc_lo < 0.0 ){
			
     		gsl_function F; 
       		F.function = &get_Sp_foc;
       		F.params   = &params;
       		
     		const gsl_root_fsolver_type *T;
	    	gsl_root_fsolver *s;
			
       		T = gsl_root_fsolver_brent;
       		s = gsl_root_fsolver_alloc (T);
       		gsl_root_fsolver_set (s, &F, x_lo, x_hi);
			
 			iter = 0; 
			max_iter = 100; 

			double eps_rel = 1.0e-4; 
 			double eps_abs  = 10.0;
//			if ( x_hi - x_lo < 100 ) eps_abs = 10.0; // 0.01 * (x_hi - x_lo); 
			
			double residual = foc_lo; 
			
 	   		do
         	{
           		iter++;
           		status = gsl_root_fsolver_iterate (s);
          		x_r  = gsl_root_fsolver_root (s);
           		x_lo = gsl_root_fsolver_x_lower (s);
           		x_hi = gsl_root_fsolver_x_upper (s);
           		
           		residual = get_Sp_foc(x_r, &params); 
           		
         		status0 = gsl_root_test_interval (x_lo, x_hi, eps_abs, eps_rel);
		   		
		   		status = gsl_root_test_residual (residual, eps_rel); 
		   		
		   		#ifdef DEBUG_Sp_foc
           		if (status == GSL_SUCCESS){ 
					printf ("\t Converged: v_r %f\n", get_value(x_r, params.cashathand-x_r, params.de, params.dk, params.dq, params.infoset));
				}   
				printf ("\t %5d [%.7f, %.7f] c_r %f x_r %.7f %.7f \t residual %f \t boundary [%.4f, %.4f]\n", iter, x_lo, x_hi, params.cashathand-x_r, x_r, x_hi - x_lo, residual, x0_lo, x0_hi);
		   		#endif
		   	
         	}
       		while (status == GSL_CONTINUE && status0 == GSL_CONTINUE && iter < max_iter);
     
    		gsl_root_fsolver_free (s);
		
			if ( status != GSL_SUCCESS && status0 != GSL_SUCCESS ) 	printf ("\t %5d [%.7f, %.7f] c_r %f x_r %.7f %.7f \t residual %f \t boundary [%.4f, %.4f]\n", iter, x_lo, x_hi, params.cashathand-x_r, x_r, x_hi - x_lo, residual, x0_lo, x0_hi);
        
        
            v_best = get_valueEGM(x_r, params.cashathand-x_r, params.de, params.dk, params.dq, params.infoset);
       
        
            if (v_best >= v0_hi -1.0e-10 && v_best >= v0_lo -1.0e-10){
                
                if ( status == GSL_SUCCESS && x_hi - x_lo < eps_abs ){   // x_lo > x0_lo + 1.0e-4 && x_hi < x0_hi - 1.0e-4
                   isbad = 0;
                }
            
            }
        
	}
	else if (foc_lo >= 0.0 && foc_hi > 0.0){ // credit constrained 
		
		x_r = x_lo;
        v_best = v0_lo; // get_value(x_r, params.cashathand-x_r, params.de, params.dk, params.dq, params.infoset);
    
       
        double x1 = (x_hi+x_lo)/2.0;
        double v1 = get_valueEGM(x1, params.cashathand-x1, params.de, params.dk, params.dq, params.infoset);
        
        if (v_best >= v0_hi && v_best >= v1){
            
            isbad = 1;
            
        } else {
            
            v_best = max( v0_hi, v1);
            
        }
    
        
	} else if (foc_lo < 0.0 && foc_hi <= 0.0){ // saving constrained
		
		x_r = x_hi;
        v_best = v0_hi; // get_value(x_r, params.cashathand-x_r, params.de, params.dk, params.dq, params.infoset);
    
       
        double x1 = (x_hi+x_lo)/2.0;
        double v1 = get_valueEGM(x1, params.cashathand-x1, params.de, params.dk, params.dq, params.infoset);
        
        if (v_best >= v0_lo && v_best >= v1){
         
            isbad = 2;
        
        } else {
            
            v_best = max( v0_lo, v1);
        }
    
    
        
	} else { // (foc_hi < 0.0 && foc_lo > 0.0 )
		
		// foc_hi < 0 while foc_lo > 0: possible due to the existence of consumption subsidy from parents. 
		// ... must make sure that x_r is between [x_lo, x_hi]
		x_r = (x_hi+x_lo)/2.0;
        v_best = get_valueEGM(x_r, params.cashathand-x_r, params.de, params.dk, params.dq, params.infoset);
				
	}
	

	
		
	if (v0_hi > v_best){
		x_r    = x0_hi; 
		v_best = v0_hi; 
	}
	if (v0_lo > v_best){
		x_r    = x0_lo; 
		v_best = v0_lo; 	
	}	

	foc_xr = get_Sp_foc(x_r, &params);  
		
	return x_r; 	
	
}
*/

void optPolicyEGM(double &max_v, int &opt_ichoice, int &opt_de, double &opt_dk, int &opt_dq, double &opt_c, double &opt_snext, double &opt_wage, double &opt_gtransfer, double &opt_gtrcmin, double &opt_ptransfer, double &opt_debtlimit_p, const struct_infoset &infoset){


	const int ia      = infoset.ia; 
	const int educ    = infoset.educ; 		
	const int ieduc   = educ - MIN_educ;
/*	
	const int itype   = infoset.isavingspr; 
	const int ieduccatpr = infoset.ieduccatpr; 
	double logcprice  = infoset.logcprice; 
*/	
	#ifdef DEBUG
	//if (kstock + ieduc > ia + (MAX_initeduc - MIN_educ) ) 
	if (infoset.kstock > ia || ieduc > ia + (MAX_initeduc - MIN_educ)) 
	{
		printf("\n\n Errors in optPolicy(): Violation in kstock %2f + ieduc %2d <= ia %2d \n ... exiting... ", infoset.kstock, ieduc, ia); 
		exit(1);  
	}
	int de_prev0 = infoset.de_prev;
	if (de_prev0 == 1 && ia-1 >= NUM_ageEduc ) {
		printf("\n\n Errors in optPolicy(): Violation in de_prev %d, ia %2d  \n ... exiting... ", de_prev0, ia); 
		exit(1);  
	}
	#endif
	
	
	//--- initialization 
	max_v        = -1.0e+30;
		
	//--- choices
	opt_ichoice  = -9; 
	opt_de       = -9; 
	opt_dk       = -9; 
	opt_dq       = -9; 
	opt_c        	= -9.0; 
	opt_snext       = -9.0; 
	
	
	//--- outcomes under choices
	opt_wage      = 0.0; 
	opt_gtransfer = 0.0; 
	opt_gtrcmin   = 0.0; 
	opt_ptransfer = 0.0; 
	opt_debtlimit_p = 0.0; 
	
	
	
	int educ_available = (ia < NUM_ageEduc ) * (ieduc < NUM_ieduc-1 ); 
	if (ia+INIT_age > MAX_ageHS && ieduc+MIN_educ < 12) educ_available = 0; 
	
	int nchoices = 2 + 4 * (educ_available==1); 
	
		
	for (int ichoice = 0; ichoice < nchoices; ichoice++){ 
		
		int de, dq; double dk; 
		revert_ichoice(ichoice, de, dq, dk); 	
		
		
		double v = -1.0e+30; double c = 0.0; double snext = 0.0; 
				
		double m, v_app, c_app, snext_app, dv, gtransfer, gtrcmin, debtlimit_p, sp_lo, consumption_min; 
		Interp_ValuePolicy_EGM_new(m, v_app, c_app, snext_app, dv, gtransfer, gtrcmin, debtlimit_p, sp_lo, consumption_min, de, dk, dq, infoset); 
		
		
		v = v_app;  c = c_app; snext = snext_app; // initialization. 
		
	if ( c > 1.0e-4 ) {  // calculate value function for those with enough resources to choose ichoice... 
		
		v = get_valueEGM(snext, c, de, dk, dq, infoset); 
	
	}
				
		
		if (max_v <= v){
				
				max_v = v; 
				
				opt_ichoice  = ichoice; 
				opt_de       = de; 
				opt_dk       = dk; 
				
				opt_dq	     = dq;
				
				opt_c        = c;
				
				opt_snext    = snext; 
				
				opt_gtrcmin   = gtrcmin; 
				opt_gtransfer = gtransfer; 
				
				opt_debtlimit_p = debtlimit_p; 
			}					
		
	} // ichoice
	
	
	
	if (opt_dk>0.1) {
		opt_wage    = get_wagerate(opt_dk, opt_de, infoset); 
	}
	opt_ptransfer  = get_ptransfer(opt_de, opt_dk, infoset); 
	


	
	#ifdef DEBUG
			
	double cost_dq = infoset.health; 

	if (opt_snext>0 && opt_c < cu_min  - cost_dq * opt_dq - 1.0e-4) printf("\n Errors in optPolicy: opt_snext = %f >0  & opt_c = %f < cu_min %f \t de %d, dk %f, dq %d \n exiting... \n", opt_snext, opt_c, cu_min, opt_de, opt_dk, opt_dq);
	
	if (max_v < 0.0) {


		printf("\n Errors in optPolicy: max_v %f \n exiting... \n", max_v); 	

		printf("\n\n ia %2d, ieduccatpr %d, isavingspr %d, de_prev %d, ieduc %2d, kstock %2f, savings %f \n", ia, infoset.ieduccatpr, infoset.isavingspr, infoset.de_prev, ieduc, infoset.kstock, infoset.savings); 
	
		printf("\n max_v %f, opt_ichoice %d, opt_de %d, opt_dk %f, opt_dq %d, opt_c %f, opt_snext %f, opt_wage %f, ia %d, ieduc %2d, kstock %2f, s %f \n", max_v, opt_ichoice, opt_de, opt_dk, opt_dq, opt_c, opt_snext, opt_wage, ia, ieduc, infoset.kstock, infoset.savings); 
		
			if (opt_c < 1000.0-10.0) {printf("\n Errors in optPolicyEGM() exiting.  c %f", opt_c); exit(1); }  // only check when iEGM == 0 


		double util = get_utilityEGM(opt_de, opt_dk, opt_c, opt_dq, ia, infoset.de_prev, ieduc, infoset.ieduccatpr, infoset.theta, infoset.qstock, infoset.eps_ue, infoset.eps_uq);
				
		printf("\n\n Errors in optPolicy(%2d): Vq = %.8f < 0.0 \t, snext %.4f, c = %.4f, choices [%d, %.1f, %d], util %f \n ... exiting ...", infoset.ia, max_v, opt_snext, opt_c, opt_de, opt_dk, opt_dq, util); 
	
		exit(1);  		
	}	
	
	#endif	
	
}


void PolicySolutionUpper(int ichoice, int ia, int de_prev, int ieduc, int itype, int ieduccatpr, int ipt){
		
//    if (diff_m < 0.0 || diff_v < 0.0){ 

		double M[NUM_MGrid], C[NUM_MGrid], Vt[NUM_MGrid]; 
    	double uM[NUM_MGrid], uC[NUM_MGrid], uVt[NUM_MGrid]; 

	double m_bigger = 1.0e+30; 	double diff_m = 10000.0; 
	double v_bigger = 1.0e+30;  double diff_v = 10000.0; 

	
    if (ia >= NUM_ageEduc) { // 
    
		int ia_ageOld = ia - NUM_ageEduc;     
		
		for (int iEGM = NUM_MGrid-1; iEGM >= 0; iEGM--){
			

			if (SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] < 0.0 || SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] > 1.0e+10
			|| SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] < 1000.0/1.2) {
			   	printf("\n Errors in PolicySolutionUpper() ia %3d, iEGM %d, M [%f, %f, %f ], c [%f, %f, %f ]", ia, iEGM, SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], M[iEGM], uM[iEGM], SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], C[iEGM], uC[iEGM]);
			}
			
		   	M[iEGM] = SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM]; 
 		   	C[iEGM] = SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM]; 
  		  	Vt[iEGM] = SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM]; 
    	
 		   	uM[iEGM] = SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM]; 
	    	uC[iEGM] = SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM]; 
 		   	uVt[iEGM] = SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM]; 
			
			// check to see if mgrid is monotonic 
			double diff_new = m_bigger - M[iEGM]; 	
			if (diff_new < diff_m) diff_m = diff_new; 
			m_bigger = M[iEGM]; 
		
			double diff_v_new = v_bigger - Vt[iEGM]; 	
			if (diff_v_new < diff_v) diff_v = diff_v_new; 
			v_bigger = Vt[iEGM]; 

			
		} // iEGM 
   		
   		if (diff_m < 0.0 || diff_v < 0.0){
   		

		UpperEnvelope(&uM[0], &uC[0], &uVt[0], &M[0], &C[0], &Vt[0], NUM_MGrid);  // specified in mexUpperEnvelope.c 
    	
    	for (int iEGM = 0; iEGM < NUM_MGrid; iEGM ++){ 
    	
			if (uM[iEGM] < 0.0 || uM[iEGM] > 1.0e+10
			|| uC[iEGM] < 1000.0/1.2) {
			   	printf("\n Errors in PolicySolutionUpper() ia %3d, iEGM %d, M [%f, %f, %f ], c [%f, %f, %f ], v [%f, %f, %f ]", ia, iEGM, SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], M[iEGM], uM[iEGM], SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], C[iEGM], uC[iEGM], SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], Vt[iEGM], uVt[iEGM]);
			}
			    	
    	SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = uM[iEGM]; 
    	SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = uC[iEGM]; 
    	SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = uVt[iEGM]; 
		
			if (SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] < 0.0 || SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] > 1.0e+10
			|| SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] < 1000.0/1.2) {
			   	printf("\n Errors in PolicySolutionUpper() ia %3d, iEGM %d, M [%f, %f, %f ], c [%f, %f, %f ]", ia, iEGM, SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], M[iEGM], uM[iEGM], SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM], C[iEGM], uC[iEGM]);
			}
			
    	} // iEGM

		} // if (diff_m < 0.0 || diff_v < 0.0)

    } // if ia >= NUM_ageEduc
    else {

		int ia_ageYoung = ia; 
			    				
    	for (int iEGM = NUM_MGrid-1; iEGM >= 0; iEGM--){


			if (SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] < 0.0 || SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] > 1.0e+10
			|| SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] < 1000.0/1.2 ) {
			   	printf("\n Errors in PolicySolutionUpper() ia %3d, iEGM %d, M_next %f, M [%f, %f, %f ], c [%f, %f, %f ], v [%f, %f, %f ]", ia, iEGM, SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM+1], SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], M[iEGM], uM[iEGM], SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], C[iEGM], uC[iEGM], SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], Vt[iEGM], uVt[iEGM]);
			}
			
		   	M[iEGM] = SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM]; 
 		   	C[iEGM] = SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM]; 
  		  	Vt[iEGM] = SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM]; 
    	
 		   	uM[iEGM] = SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM]; 
	    	uC[iEGM] = SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM]; 
 		   	uVt[iEGM] = SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM]; 
			
			
			// check to see if mgrid is monotonic 
			double diff_new = m_bigger - M[iEGM]; 	
			if (diff_new < diff_m) diff_m = diff_new; 
			m_bigger = M[iEGM]; 
		
			double diff_v_new = v_bigger - Vt[iEGM]; 	
			if (diff_v_new < diff_v) diff_v = diff_v_new; 
			v_bigger = Vt[iEGM]; 

		} // iEGM 
   		

		
		if (diff_m < 0.0 || diff_v < 0.0){
		   		
		UpperEnvelope(&uM[0], &uC[0], &uVt[0], &M[0], &C[0], &Vt[0], NUM_MGrid);  // specified in mexUpperEnvelope.c 
    	
    	for (int iEGM = 0; iEGM < NUM_MGrid; iEGM ++){ 
			
	
			if (uM[iEGM] < 0.0 || uM[iEGM] > 1.0e+10
			|| uC[iEGM] < 1000.0/1.2) {
			   	printf("\n Errors in PolicySolutionUpper() ia %3d, iEGM %d, M [%f, %f, %f ], c [%f, %f, %f ]", ia, iEGM, SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], M[iEGM], uM[iEGM], SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], C[iEGM], uC[iEGM]);
			}
			
			    	
    	SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = uM[iEGM]; 
    	SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = uC[iEGM]; 
    	SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = uVt[iEGM]; 

	
			if (SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] < 0.0 || SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] > 1.0e+10
			|| SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] < 1000.0/1.2) {
			   	printf("\n Errors in PolicySolutionUpper() ia %3d, iEGM %d, M [%f, %f, %f ], c [%f, %f, %f ]", ia, iEGM, SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], M[iEGM], uM[iEGM], SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM], C[iEGM], uC[iEGM]);
			}

    	} // iEGM
		
		} // if (diff_m < 0.0 || diff_v < 0.0)
				
    } // if ia < NUM_ageEduc
   		 


}



void PolicySolutionEGM(int ichoice, int ia, int de_prev, int ieduc, int itype, int ieduccatpr, int ipt){
	
	
	int de; int dq; double dk; 
	revert_ichoice(ichoice, de, dq, dk); 
	
	
	double state[DIM_smolyak]; 								
	for (int idim=0; idim<DIM_smolyak; idim++){		
	state[idim] = (GRID_smolyak[ipt][idim] + 1.0)*(MAX_coefEmaxstate[ia][ieduc][idim]-MIN_coefEmaxstate[ia][ieduc][idim])/2.0 + MIN_coefEmaxstate[ia][ieduc][idim];
	}
				
	double theta[DIM_theta]; double logcprice, qstock; 				
	get_stateMODEL(state, theta[0], theta[1], logcprice, qstock, ia, ieduc); 

	    
	const int educ    = ieduc + MIN_educ; 
	double kstock = get_kstock(educ, ia); 
		
	const int 	 ieduc_next  = ieduc + de; 
	const double kstock_next = get_knext(ia, educ, kstock, dk, de);	
	double qstock_next = get_qnext_itype(qstock, dq, ia, itype);

	
    double debtlimit_p = get_debtlimit_ptotal(theta, ia, educ, de, dq, qstock); 
	double sp_lo = get_minS_EGM(ia, educ, de, dk, debtlimit_p); 
//	m = get_totalwealthEGM(de, dk, dq, infoset, gtrcmin) - sp_lo; 
	
    
	double c_ptsubsidy = getc_subsidy(ia, educ, de); 
	double c_minexpenditure = getc_minexpenditure(ia,  educ, de, c_ptsubsidy);
	double consumption_min  = max(c_ptsubsidy, c_minexpenditure); 
    
    double probs = get_survivalProb(ia, educ, qstock, dq, de); 
	double delta = get_delta(theta, qstock); // infoset.delta; 
	
	double factor_deltaprobs = pow( probs * delta, -1.0/cu_g ); 

    
	double evnext_smin; 
	
	double eps_ue = 0.0; double eps_uq = 0.0; 
	
	double m_bigger = 1.0e+30; 	double diff_m = 10000.0; 
	double v_bigger = 1.0e+30;  double diff_v = 10000.0; 
	
	
	double du_highS = 0; 
	
for (int iEGM = NUM_MGrid-1; iEGM >= NumVfi; iEGM--){
			
	double snext = Snext_plusGrid[ia][ieduc][iEGM-NumVfi] + sp_lo;  // Splus = snext - sp_lo >= 0  
		
	double loghnext_expected = logcprice; 
	{ 	
		loghnext_expected = get_logcprice_next(logcprice, 0.0); // infoset.eps_h
	} 
	

	
	// Calculate the expected marginal utility and expected value function	
	double EdV_next  = -9.0;  	double ev_next = -9.0; 
	
	if (ia < NUM_age-1){ 
		
	ev_next = get_EvnextEGM_new(theta, ia+1, ieduccatpr, itype, de, ieduc_next, kstock_next, snext, loghnext_expected, qstock_next, EdV_next); 	 
		
	} else {
	
	ev_next  = get_PDV_EGM_new(ia+1, ieduc_next, theta, kstock_next, exp(loghnext_expected), snext, itype, dq, qstock_next, EdV_next); 
	
	}
	
	if (EdV_next < du_highS ) EdV_next = du_highS; 
	


    double Inv_dVnext = inv_dUdC_EGM( EdV_next, ia, educ);  // a function of current M at the optimal 
	double c = factor_deltaprobs * Inv_dVnext; 
	if (c < consumption_min ) c = consumption_min; 
	
	
	double m = (snext - sp_lo) + c;                 						// Splus_Grid >=0, c >= c_minexpenditure, hence m >= c_minexpenditure. 
	
	if (c < 1000.0) {printf("\n Errors: c %f, consumption_min %f ", c, consumption_min); exit(1); }  // only check when iEGM == 0 
	
	double v = get_utilityEGM(de, dk, c, dq, ia, de_prev, ieduc, ieduccatpr, theta, qstock, eps_ue, eps_uq) + probs * delta * ev_next; 
	
	if (ia >= NUM_ageEduc) {
	
	int ia_ageOld = ia - NUM_ageEduc; 
	
	SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = m; 	
	SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = c/factor_deltaprobs; 
	SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = v; 
	

	} else {
	
	int ia_ageYoung = ia; 
	
	SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = m; 	
	SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = c/factor_deltaprobs; 
	SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = v; 
	
	}
	
	du_highS = EdV_next;
			
	if (iEGM == NumVfi) evnext_smin = ev_next; 
	
	// check to see if mgrid is monotonic 
	double diff_new = m_bigger - m; 	
	if (diff_new < diff_m) diff_m = diff_new; 
	m_bigger = m; 
	
	
	double diff_v_new = v_bigger - v; 	
	if (diff_v_new < diff_v) diff_v = diff_v_new; 
	v_bigger = v; 

} // iEGM, iEGM >= NumVfi
	
	// the endougenous current period cash when snext = - deblimit_sp under consumption smoothing (if current consumption > c_subsidy)
	const double m_lo = m_bigger; // SolutionM_EGM[ichoice][ia][de_prev][ieduc][itype][ieduccatpr][ipt][NumVfi]; 	
    if (m_lo < consumption_min ){ printf("\n Errors in PolicySolutionEGM(): m_lo %f < c_min %f \n Exiting ...\n", m_lo ,consumption_min); exit(1);  }   
	
    // Value function and lower part of resources (where a<=0), not found by EGM
    // Add grid of resources below M_lo, that is not solved by EGM. This is
    // done to get accurate value function. 
	
	
for (int iEGM = NumVfi-1; iEGM >= 0; iEGM--){
		
		double m = 1.0e-4; 
		double c = m;
		{
			m = consumption_min + 1.0*iEGM/NumVfi * (m_lo - consumption_min);   // m > c_minexpenditure
			c = m; 
		} 
		
		
	if (c < 1000.0) {printf("\n Errors!!! exiting.  c %f, consumption_min %f ", c, consumption_min); exit(1); }  // only check when iEGM == 0 
	
	double v = get_utilityEGM(de, dk, c, dq, ia, de_prev, ieduc, ieduccatpr, theta, qstock, eps_ue, eps_uq) + probs * delta * evnext_smin; 
	
	
	if (ia >= NUM_ageEduc) {

	int ia_ageOld = ia - NUM_ageEduc; 
	
	SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = m; 	
	SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = c/factor_deltaprobs; 
	SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][iEGM] = v; 

		
	} else {
	
	int ia_ageYoung = ia; 
	
	SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = m; 	
	SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = c/factor_deltaprobs; // c; 
	SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][iEGM] = v; 
	
	}	
				
		
		// check to see if mgrid is monotonic 
		double diff_new = m_bigger - m; 	
		if (diff_new < diff_m) diff_m = diff_new; 
		m_bigger = m; 
		
		double diff_v_new = v_bigger - v; 	
		if (diff_v_new < diff_v) diff_v = diff_v_new; 
		v_bigger = v; 
		
} // iEGM, iEGM < NumVfi


	#ifdef DEBGU	
	if ( (diff_m < 0 || diff_v < 0 ) && ia == NUM_age-1) {	printf("\n Error: last period cash at hand is not monotonic: diff_m %f , diff_v %f  \n", diff_m, diff_v ); exit(1); }
	#endif 

	
	
}











//  Solving the model
void SolveModel_EGM(){

			const int minitype = 0; 
			const int maxitype = 0; 
		
     
    gsl_vector * coefValue = gsl_vector_alloc(NUM_coefEmax); 
	struct_smolyak smolyak; 

	       
	       
            //  Initialize struct to hold NumChoice choices
            for (int ichoice=0; ichoice < NumChoice_ageYoung; ichoice++){
            
            //!!!! Initialization of value function coefficients
			for (int ia = 0; ia < NUM_ageEduc; ia++){
			
			for (int de_prev = 0; de_prev < NUM_de; de_prev++){  // whether left school or not
			for (int ieduc = 0; ieduc < NUM_ieduc; ieduc++){

			for (int itype = minitype; itype <= maxitype; itype++){
			
			
			for (int ieduccatpr = 0; ieduccatpr < NUM_ieduccatpr; ieduccatpr++){
			
			
			for (int im = 0; im < NUM_MGrid; im++){
				
				for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
				
				// FOR each ipt in sparse grid, calculate endogenous M, C, and V.
				SolutionM_ageYoung_EGM[ichoice][ia][de_prev][ieduc][itype][ieduccatpr][ipt][im] = 0.0; 
				SolutionC_ageYoung_EGM[ichoice][ia][de_prev][ieduc][itype][ieduccatpr][ipt][im] = 0.0; 
				SolutionV_ageYoung_EGM[ichoice][ia][de_prev][ieduc][itype][ieduccatpr][ipt][im] = 0.0; 
								
				// For each m in Common_MGrid, interpolation Coef of sparse grid [thetac, thetan, logh, q]				
//				CoefC_ageYoung_EGM[ichoice][ia][de_prev][ieduc][itype][ieduccatpr][im][ipt] = 0.0; 
				CoefV_ageYoung_EGM[ichoice][ia][de_prev][ieduc][itype][ieduccatpr][im][ipt] = 0.0; 
				
				} // ipt 
							
			} // grid for m 
			
			}}
			}}
			}
			} // ichoice
			
			
			//  Initialize struct to hold NumChoice choices
            for (int ichoice=0; ichoice < NumChoice_ageOld; ichoice++){
            
            //!!!! Initialization of value function coefficients
			for (int ia = 0; ia < NUM_ageNoEduc; ia++){
			
//			for (int de_prev = 0; de_prev < NUM_de; de_prev++){  // whether left school or not
			for (int ieduc = 0; ieduc < NUM_ieduc; ieduc++){
						
			for (int itype = minitype; itype <= maxitype; itype++){
						
			
			for (int im = 0; im < NUM_MGrid; im++){
				
				for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
				
				// FOR each ipt in sparse grid, calculate endogenous M, C, and V.
				SolutionM_ageOld_EGM[ichoice][ia][ieduc][itype][ipt][im] = 0.0; 
				SolutionC_ageOld_EGM[ichoice][ia][ieduc][itype][ipt][im] = 0.0; 
				SolutionV_ageOld_EGM[ichoice][ia][ieduc][itype][ipt][im] = 0.0; 
				
				
				// For each m in Common_MGrid, interpolation Coef of sparse grid [thetac, thetan, logh, q]				
//				CoefC_ageOld_EGM[ichoice][ia][ieduc][itype][im][ipt] = 0.0; 
				CoefV_ageOld_EGM[ichoice][ia][ieduc][itype][im][ipt] = 0.0; 
								
				} // ipt 
							
			} // grid for m 
			
			}}
//			}}
			}
			} // ichoice
			
			
            
		for (int ia = NUM_age-1; ia >= NUM_ageEduc; ia--){
			
			int ia_ageOld = ia - NUM_ageEduc; 
			//--- only calculate for possible schooling and experience levels for given age				
			int maxieduc  = min( NUM_ieduc-1, ia + (MAX_initeduc - MIN_educ) ); 	
			
			
			for (int ieduc = 0; ieduc <= maxieduc; ieduc++){ 
			
			for (int itype = minitype; itype <= maxitype; itype++){
			
			
			int de_prev = 0; int ieduccatpr = 0; 
			

			#ifdef DEBUG_coefEmax
	    	printf("\n emax ia %3d, age %d, ieduc %2d, itype %d \t", ia, ia+INIT_age, ieduc, itype); 
	    	#endif
			
					
			for (int ichoice = 0; ichoice < NumChoice_ageOld; ichoice++){ 
			
			
				
				for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
					
					// for each ipt, update C, V, M for a vector of cash-at-hand NUM_MGrid. 
					PolicySolutionEGM(ichoice, ia, de_prev, ieduc, itype, ieduccatpr, ipt); 
					
				} // ipt 
								
				#ifdef USE_omp
			
				#pragma omp barrier
				#pragma omp flush
			
				#endif
				
				

				// check the optimal solution 
//				#ifdef MODEL_upper
				


				for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
					
					
					PolicySolutionUpper(ichoice, ia, de_prev, ieduc, itype, ieduccatpr, ipt); 
					
						for (int im = 1; im < NUM_MGrid; im++){  
							double diff_m_new = SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] - SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im-1]; 
							double diff_v_new = SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] - SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im-1]; 
							// SolutionV may be negative due to utility(ichoice)
							if (SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] < 0.0 || SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] < 0.0 || SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] < 0.0 || diff_m_new < 0.0 || diff_v_new < 0.0 || SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] < 0.0 || SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im] < 0.0 ) {
								printf("\n Errors in SolveModel_EGM(): itype %d, ia %d, ieduc %d, ipt %d, im %d, m = %f, m_prev %f, v = %f, v_prev = %f, c=%f, c_prev = %f, \n exiting \n", itype, ia, ieduc, ipt, im, SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im], SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im-1], SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im], SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im-1], SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im], SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im-1]); 
								exit(1); 
							}
						}	

					
					
				} // ipt
				
//				#endif 
													
				
																
				#ifdef USE_omp
			
				#pragma omp barrier
				#pragma omp flush
			
				#endif
				

				double vec_C[NUM_MGrid][NUM_coefEmax], vec_V[NUM_MGrid][NUM_coefEmax]; 

				
				// update the Smolyak coefficient at each M grid 
				for (int im = 0; im < NUM_MGrid; im++){ 
					
					// a common grid that is invariant to the sparse grid points: Common_MGrid[ichoice][ia][ieduc][itype][NUM_MGrid]
					// initialization using ipt=0					
					Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im] = SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][0][im]; 
						
	
					vec_C[im][0]                                                  = SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][0][im];
					vec_V[im][0]                                                  = SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][0][im]; 
					for (int ipt = 1; ipt < NUM_coefEmax; ipt++){
						
						// value of C at ipt by interpolating over SolutionM 
						LinInterp(&vec_C[im][ipt], &SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][0], &SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][0], &Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im], NUM_MGrid,1);
						
						// value of V at ipt by interpolating over SolutionM 
						LinInterp(&vec_V[im][ipt], &SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][0], &SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][0], &Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im], NUM_MGrid,1);
						
						
//						vec_logCMratio[im][ipt] = log ( vec_C[im][ipt] / Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im] ) ; 
						
						if ( vec_C[im][ipt] <= 1.0e-6 || Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im]  <= 1.0e-6 ) {
							printf("\n Errors in SolveModel_EGM(): c <= 0.0 or m <= 0.0 \t ichoice %d, ipt %d, im %d, common m %f, v %f, c %f, \t endogenous m %f, v %f, c %f ", ichoice, ipt, im, Common_ageOld_MGrid[ichoice][ia_ageOld][ieduc][itype][im], vec_V[im][ipt], vec_C[im][ipt], SolutionM_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im], SolutionV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im], SolutionC_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][ipt][im]); 
							exit(1); 
						}
						
					}
					
					for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
						Common_ageOld_optC[ichoice][ia_ageOld][ieduc][itype][im][ipt] = vec_C[im][ipt];
					}
					
					// V 
					smolyak.PolyCoef(ChebyIndex, POLYINV_smolyak, vec_V[im], coefValue); 	
					for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
						CoefV_ageOld_EGM[ichoice][ia_ageOld][ieduc][itype][im][ipt]  = gsl_vector_get(coefValue, ipt);	
					}
					
				} // im 
				
				
												
				#ifdef USE_omp
			
				#pragma omp barrier
				#pragma omp flush
			
				#endif

		} // ichoice
		
						
		}} // ieduc, itype
		
	} // ia >= NUM_ageEduc
	
    	
		for (int ia = NUM_ageEduc-1; ia >= 0; ia--){
			
			int ia_ageYoung = ia ; 
							
			//--- only calculate for possible schooling and experience levels for given age				
			int maxieduc  = min( NUM_ieduc-1, ia + (MAX_initeduc - MIN_educ) ); 	
								
			for (int ieduc = 0; ieduc <= maxieduc; ieduc++){ 
			
			// ... de_prev only matters for current schooling choices ... 
			int n_de   = 2;  int nchoices = 6; 
				
			int educ_available = (ia < NUM_ageEduc ) * (ieduc < NUM_ieduc-1 ); 
			if (ia+INIT_age > MAX_ageHS && ieduc+MIN_educ < 12) educ_available = 0; 
			if (educ_available==0){ n_de = 1; nchoices = 2; }
			 
			for (int de_prev=0; de_prev < n_de; de_prev++){
						
			int n_ieduccatpr  = 1 * (ia >= NUM_age_ptr) + NUM_ieduccatpr * (ia < NUM_age_ptr); 
			for (int ieduccatpr = 0; ieduccatpr < n_ieduccatpr; ieduccatpr++){
			
			
			for (int itype = minitype; itype <= maxitype; itype++){
			
			
			#ifdef DEBUG_coefEmax
	    	printf("\n emax ia %3d, de_prev %d, ieduc %2d, itype %d \t", ia, de_prev, ieduc, itype); 
	    	#endif
			
			
			for (int ichoice = 0; ichoice < nchoices; ichoice++){ 
			
			
				
				for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
					
					// for each ipt, update C, V, M for a vector of cash-at-hand NUM_MGrid. 
					PolicySolutionEGM(ichoice, ia, de_prev, ieduc, itype, ieduccatpr, ipt); 
									
				} // ipt 
								
				#ifdef USE_omp
			
				#pragma omp barrier
				#pragma omp flush
			
				#endif
				
				
				
				// check the optimal solution 
//				#ifdef MODEL_upper				
				
				
				
				for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
					
					PolicySolutionUpper(ichoice, ia, de_prev, ieduc, itype, ieduccatpr, ipt); 								

						for (int im = 1; im < NUM_MGrid; im++){  
							double diff_m_new = SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im] - SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im-1]; 
							double diff_v_new = SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im] - SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im-1]; 
							
							if (diff_m_new < 0.0 || diff_v_new < 0.0 || SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im] < 0.0 || SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im] < 0.0 ) {
								printf("\n Errors in SolveModel_EGM(): ipt %d, im %d, m [%f, %f], v [%f, %f], c [%f, %f], \n exiting \n", ipt, im, SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im-1], SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im], SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im-1], SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im], SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im-1], SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im]); 
								exit(1); 
							}
						}	

				
				} // ipt 
				
//				#endif 		
				
				#ifdef USE_omp
			
				#pragma omp barrier
				#pragma omp flush
			
				#endif
				
								
					
				// update the Smolyak coefficient at each M grid 
				for (int im = 0; im < NUM_MGrid; im++){ 
					
					// a common grid that is invariant to the sparse grid points: Common_MGrid[ichoice][ia][ieduc][itype][NUM_MGrid]
					// initialization using ipt=0					
					Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im] = SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][0][im]; 
					
										
					double vec_C[NUM_coefEmax], vec_V[NUM_coefEmax]; 
									
					vec_C[0]                                                  							= SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][0][im];
					vec_V[0]                                                  							= SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][0][im];
					for (int ipt = 1; ipt < NUM_coefEmax; ipt++){
						
						// value of C at ipt by interpolating over SolutionM 
						LinInterp(&vec_C[ipt], &SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][0], &SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][0], &Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im], NUM_MGrid,1);
						
						// value of V at ipt by interpolating over SolutionM 
						LinInterp(&vec_V[ipt], &SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][0], &SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][0], &Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im], NUM_MGrid,1);
						

						
						
						if ( vec_C[ipt] <= 1.0e-6 || Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im]  <= 1.0e-6 ) {
							printf("\n Errors in SolveModel_EGM(): c <= 0.0 or m <= 0.0 \t ichoice %d, ipt %d, im %d, common m %f, v %f, c %f, \t endogenous m %f, v %f, c %f ", ichoice, ipt, im, Common_ageYoung_MGrid[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im], vec_V[ipt], vec_C[ipt], SolutionM_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im], SolutionV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im], SolutionC_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][ipt][im]); 
							exit(1); 
						}
					
					}

					for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
						Common_ageYoung_optC[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im][ipt] = vec_C[ipt];
					}

			
					// V 
					smolyak.PolyCoef(ChebyIndex, POLYINV_smolyak, vec_V, coefValue); 	
					for (int ipt = 0; ipt < NUM_coefEmax; ipt++){
						CoefV_ageYoung_EGM[ichoice][ia_ageYoung][de_prev][ieduc][itype][ieduccatpr][im][ipt]  = gsl_vector_get(coefValue, ipt);	
					}
					
					
					
				} // im 
								
				#ifdef USE_omp
			
				#pragma omp barrier
				#pragma omp flush
			
				#endif
				
		} // ichoice 
				
		}}}}
         
        }   
    

}

	
void get_smOpt_EGM(int miniper, int maxiper, int minia, int maxia){
	
   
	for (int ia = minia; ia <= maxia; ia++){ // NUM_age
	
	
	
	
	for (int iper = miniper; iper <= maxiper; iper++){	
		
		
		double theta[] = {sm_theta[0][iper], sm_theta[1][iper]};
		
			sm_isavingspr[iper]  =  0; 

			sm_wgt[iper]          = 1.0; 
					
		struct_infoset infoset;
		
		infoset.theta = theta; 

		infoset.health   = sm_cprice[ia][iper]; 
		infoset.logcprice   = log(sm_cprice[ia][iper]); 
				
 		infoset.qstock   = sm_qstock[ia][iper];  
		infoset.delta   = get_delta(infoset.theta, infoset.qstock); 						
				
		
		
		infoset.ia       = ia; 
		infoset.de_prev  = sm_de_prev[ia][iper];    // left school at age ia-1
		
		infoset.debtlimitLag = 0.0; 
		if (ia>0) infoset.debtlimitLag = sm_debtlimit_p[ia-1][iper];
		
		infoset.educ     = sm_educ[ia][iper]; 
		infoset.kstock   = sm_kstock[ia][iper]; 

								    
		infoset.savings  = sm_savings[ia][iper]; 
		
		infoset.ieduccatpr  = sm_ieduccatpr[iper]; 		
		infoset.isavingspr  = sm_isavingspr[iper];  
		
		infoset.eps_ue 	        = sm_rng_eps[ia][iper][1] ; 
		infoset.eps_uq          = sm_rng_eps[ia][iper][2] ; 
		
		
		infoset.eps_logwf       = get_rng_epsw(infoset.educ, sm_rng_eps[ia][iper][6]); 
		
		infoset.eps_health_next = get_eps_health_next(log(sm_cprice[ia][iper]), log(sm_cprice[ia+1][iper])); // sm_rng_eps[ia][iper][7];  
		
		#ifndef INFOSET_hnext
						
		infoset.eps_health_next = 0.0; 
						
		#endif 
		
				
		infoset.eps_ptr   		= sm_rng_eps[ia][iper][8]; 	 
		infoset.eps_ptrdummy   	= sm_rng_eps[ia][iper][9]; 
		
		infoset.skill           = get_skill(sm_educ[ia][iper], theta, ia, infoset.eps_logwf); 
		
		double maxV = -1.0e+30; int opt_ichoice, opt_de, opt_dq; double opt_dk, opt_c, opt_snext, opt_wage, opt_gtransfer, opt_gtrcmin, opt_ptransfer, opt_debtlimit_p;  
			
		optPolicyEGM(maxV,opt_ichoice, opt_de, opt_dk, opt_dq, opt_c, opt_snext, opt_wage, opt_gtransfer, opt_gtrcmin, opt_ptransfer, opt_debtlimit_p, infoset); 
		
		
		
		sm_value[ia][iper]   = maxV; 		
//		sm_ichoice[ia][iper] = opt_ichoice; 
		sm_de[ia][iper]      = opt_de; 
		sm_dk[ia][iper]      = opt_dk; 
		sm_dq[ia][iper]      = opt_dq; 
		sm_c[ia][iper]       = opt_c;
	
		sm_gtransfer[ia][iper] = opt_gtransfer; 
		sm_gtrcmin[ia][iper]   = opt_gtrcmin; 
		sm_ptransfer[ia][iper] = opt_ptransfer; 		

		double avg_epsw = get_avg_epsw(sm_educ[ia][iper]); 
		sm_skill [ia][iper]  = get_skill(sm_educ[ia][iper], theta, ia, avg_epsw);  
		
		sm_work_part[ia][iper]     = (sm_dk[ia][iper]  > 0.1) * (sm_dk[ia][iper]  < 0.6); // dk=0.5
		sm_work_full[ia][iper]     = (sm_dk[ia][iper]  > 0.6); 		                      // dk=1
		sm_work_fullpart[ia][iper] = (sm_dk[ia][iper]  > 0.1); 
		
		
		sm_wage[ia][iper]      = opt_wage; 
		
		double dk = sm_dk[ia][iper]; 
        double hrs = 0.0 + HOURSP * (dk>0.1 && dk < 0.9) + HOURS * (dk > 0.9);
        sm_earnings[ia][iper]  = opt_wage * hrs;		// true earnings, free of measurement errors
        sm_hours[ia][iper]     = hrs;
        
		// update individual's state variables: 
		sm_de_prev[ia+1][iper] = sm_de[ia][iper]; 
		sm_dq_prev[ia+1][iper] = sm_dq[ia][iper]; 
		
		sm_educ [ia+1][iper]   = sm_educ [ia][iper] + sm_de[ia][iper]; 		
		
		sm_kstock [ia+1][iper] = get_knext(ia, sm_educ[ia][iper], sm_kstock [ia][iper], sm_dk[ia][iper], sm_de[ia][iper]); 
		sm_expr [ia+1][iper]   = sm_expr [ia][iper] + sm_work_fullpart[ia][iper] * (sm_de[ia][iper]==0);	
			
		sm_add  [ia+1][iper]   = sm_add  [ia][iper] + sm_dq[ia][iper]; 	
		sm_add0 [ia+1][iper]   = (sm_add[ia+1][iper] <= 0.01); 
				
		sm_qstock[ia+1][iper]  = get_qnext_itype(sm_qstock[ia][iper], sm_dq[ia][iper], ia, -1);


		sm_savings[ia+1][iper]= opt_snext;  
		

		sm_debtlimit_p[ia][iper] = opt_debtlimit_p;  
		
		sm_epslogw[ia][iper] = infoset.eps_logwf/avg_epsw; 
		
		
		sm_probs[ia][iper] = 	get_survivalProb(ia, sm_educ[ia][iper], sm_qstock[ia][iper], sm_dq[ia][iper], sm_de[ia][iper]); 
		
		
		#ifdef DEBUG_smIperOpt
		double add = 0.0; double health = 1.0; double dq = 0; 
			printf("\n\n ia %2d iper %3d, maxV %.4f, qstock %f, educ %2d, kstock %2f, add %2f, health %f, savings %f, \t opt_ichoice %d, opt_dq %d, opt_snext %.4f, opt_wage %.4f, skill %.4f \n", ia, iper, maxV, sm_qstock[ia][iper], sm_educ[ia][iper], sm_expr[ia][iper], add, health, sm_savings[ia][iper], opt_ichoice, dq, opt_snext, opt_wage, sm_skill [ia][iper]); 
		#endif 
		
		
	// double check the following condition in optPolicy().
	#ifdef DEBUG
		if (sm_value[ia][iper] < 0.0 ) {
			printf("\n\n Errors in sm_value[%2d][%6d]): Vq = %.8f < 0.0  \n ... exiting ...", ia, iper, sm_value[ia][iper]); exit(1);
		}	
	#endif 		
		
		
					
	} // iper	
	
		
		#ifdef USE_omp
		
		#pragma omp barrier
		#pragma omp flush
		
		#endif

			
	}//ia
	
	
	//!!!!!!!! Individual optimization ends here !!!!!!!!!//
}



double smm_obj_est(const gsl_vector *x, void *void_params){
	
	#ifdef USE_omp
			
	#pragma omp barrier
	#pragma omp flush
	
	#endif
		
	#ifdef DEBUG
	for (int i_s =0; i_s<Nvar; i_s++) printf("\n \t x [%3d] = %5.6f ; ", i_s, gsl_vector_get (x, i_s));
	#endif 

	getpar(x); // get the structure parameters from the input vector


	#ifdef DEBUG_coefEmax
	time_t start0,end0;	double timediff; 
	time (&start0);	
	#endif 
	
	SolveModel_EGM();

	
	
	#ifdef DEBUG_coefEmax
	time (&end0);	timediff = difftime (end0,start0);
	printf ("\n\n\t It took %f minutes to get_coefEmax() \n\n", timediff/60 );
	#endif 		
	
	
//	init_sm0(0, NUM_smmPerson-1); 	
	get_smOpt_EGM(0, NUM_smmPerson-1, 0, 35-INIT_age); 
	
	
    double sumf2 = get_sumf2();
	
	
	#ifdef USE_omp
			
	#pragma omp barrier
	#pragma omp flush
	
	#endif
	
	return sumf2; 
}

double smm_obj(const gsl_vector *x, void *void_params){
	
	#ifdef USE_omp
			
	#pragma omp barrier
	#pragma omp flush
	
	#endif
		
	#ifdef DEBUG
	for (int i_s =0; i_s<Nvar; i_s++) printf("\n \t x [%3d] = %5.6f ; ", i_s, gsl_vector_get (x, i_s));
	#endif 
	
	getpar(x); // get the structure parameters from the input vector



	#ifdef DEBUG_coefEmax
	time_t start0,end0;	double timediff; 
	time (&start0);	
	#endif 
	
	SolveModel_EGM();

	
	
	#ifdef DEBUG_coefEmax
	time (&end0);	timediff = difftime (end0,start0);
	printf ("\n\n\t It took %f minutes to get_coefEmax() \n\n", timediff/60 );
	#endif 		
	
	
//	init_sm0(0, NUM_smmPerson-1); 	
	get_smOpt_EGM(0, NUM_smmPerson-1, 0, NUM_age-1); 
	
	
    double sumf2 = get_sumf2();
	
	
	#ifdef USE_omp
			
	#pragma omp barrier
	#pragma omp flush
	
	#endif
	
	return sumf2; 
}


double get_sumf2(){
    
	double sumf2 = 0.0; 	
	double  f[NUM_ismm];  
	double  moments_vec[NUM_ismm]; 
	for (int ismm=0; ismm < NUM_ismm; ismm++){
		
		moments_vec[ismm] = 0.0; 
		
		// check code, set the default to be target mean
		moments_vec[ismm] = TARGET_moments[ismm]; 

	}

	
	#ifdef USE_omp
			
	#pragma omp barrier
	#pragma omp flush
	
	#endif
	

	int NUM_smmAge = NUM_age30; 
	
	
	int age_savings[] = {20, 25, 30}; 
	int numcat_savings = 3; 
	
	//.. education at age 25 
	int age2530[] = {25, 30}; 
    
    
    
	
			
	for (int iper =0; iper < NUM_smmPerson; iper++){
				
		for (int i2530 = 0; i2530 < 2; i2530++) { 	
			
			int ia = age2530[i2530] - INIT_age; 
			
			sm_educcat[i2530][0][iper] = (sm_educ[ia][iper] < 12 ); 
			sm_educcat[i2530][1][iper] = (sm_educ[ia][iper] == 12 ); 
			sm_educcat[i2530][2][iper] = (sm_educ[ia][iper] > 12 ) * (sm_educ[ia][iper] < 16 );
			sm_educcat[i2530][3][iper] = (sm_educ[ia][iper] >= 16 ); 
		}	
	
	}


	
	int ivec = 0; 
	
	for (int ia = 16-INIT_age; ia <= 30-INIT_age; ia++){
		
        moments_vec[ivec] = my_mean(NUM_smmPerson, sm_educ[ia]);
		
        ivec++;

	}// iper
	
	
	for (int ia = 16-INIT_age; ia <= 30-INIT_age; ia++){

		moments_vec[ivec] = my_mean(NUM_smmPerson, sm_add[ia]);  
		ivec++; 
			
	}// iper
	

	
    // ... part-time working and in school
    int n_workpart_de1 = 0; int n_workpart_de1_tot = 0;
    
    // ... part-time working and in college 
    n_workpart_de1 = 0; n_workpart_de1_tot = 0;
    double n_earnings_de1 = 0.0; int n_earnings_de1_tot = 0; 
    
    int n_preduc = 0; int n_delag = 0; 
    
	for (int ia = 0; ia <= 29-INIT_age; ia++){ 
	    for (int iper =0; iper < NUM_smmPerson; iper++){
            
            if (sm_de[ia][iper] == 1 && sm_educ[ia][iper] >= 12 && sm_educ[ia][iper] < 16 ) {
                
                n_workpart_de1     += sm_work_part[ia][iper];
                n_workpart_de1_tot += 1;
                
                n_preduc += (sm_ieduccatpr[iper]==1); 
                n_delag  += sm_de_prev[ia][iper]; 
                
                if (sm_earnings[ia][iper] > 1.0e-6){
                    n_earnings_de1     += log(sm_earnings[ia][iper]);
                    n_earnings_de1_tot += 1;
                }
                
            }
            
        }// age <= 29 
    }// iper

    moments_vec[ivec] = 1.0 * n_preduc/n_workpart_de1_tot;
    ivec++;

    moments_vec[ivec] = 1.0 * n_delag/n_workpart_de1_tot;
    ivec++;
    
    moments_vec[ivec] = 1.0 * n_workpart_de1/n_workpart_de1_tot;
    ivec++;
    
    moments_vec[ivec] = 1.0 * n_earnings_de1/n_earnings_de1_tot;
    ivec++;
	
	
	int ia25 = 25 - INIT_age; 
	
	int ndq = 33-INIT_age-MIN_smmAge+1; 
	
	double vv_dq_clg[ndq], vv_dq_nclg[ndq]; 
	int    sm_tot_add0[ndq]; int    sm_dq_add0[ndq];	
	
	for (int ia = MIN_smmAge; ia < MIN_smmAge + ndq; ia++){
		
	
		
		int tot_add0=0; int tot_add0_dq=0; 
		double n_clg = 0.0; double sm_q_clg = 0.0; double sm_q_nclg = 0.0; 
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
			tot_add0    += sm_add0[ia][iper]; 
			tot_add0_dq += sm_dq[ia][iper] * sm_add0[ia][iper]; 
			
			if (sm_educ[ia25][iper] >= 16){
				
				n_clg ++; 
				
				sm_q_clg += sm_dq[ia][iper]; 
				
			} else {
				
				sm_q_nclg += sm_dq[ia][iper]; 	
			}
		
		} // iper 
		
		sm_tot_add0[ia-MIN_smmAge] = tot_add0;  
		sm_dq_add0[ia-MIN_smmAge]  = tot_add0_dq; 
		
		vv_dq_nclg[ia-MIN_smmAge] = sm_q_nclg/(NUM_smmPerson-n_clg); 

		vv_dq_clg[ia-MIN_smmAge] = sm_q_clg/n_clg; 
				
	}
	
	
	
	// smoking initiation ratio: from age 15 to age 29 (age 30 initiation rate is not available)
	for (int ia = MIN_smmAge; ia <= 29 - INIT_age; ia++){
				
		moments_vec[ivec] = 0.0; 
		if (sm_tot_add0[ia]>=1){
		moments_vec[ivec] = 1.0*sm_dq_add0[ia]/sm_tot_add0[ia]; 
		}
		ivec++; 	
		
	}
	

	//... log earning when not in school ... 
	for (int ia = 16-INIT_age; ia < NUM_smmAge; ia++){

		double sum_logearnings = 0.0; double n_logearnings = 0.0; 
		for (int iper =0; iper < NUM_smmPerson; iper++){
		
			if (sm_de[ia][iper]==0 && sm_earnings[ia][iper] > 0.0){
			
			sum_logearnings += log(sm_earnings[ia][iper]); 
			n_logearnings   += 1.0; 
			
			}
		}
		
		if (n_logearnings < 1.0 ) n_logearnings=1.0; 
		
		moments_vec[ivec] = sum_logearnings/n_logearnings;  
		ivec++;
		
	}
			

	
	for (int i2530 = 1; i2530 < 2; i2530++) { 	
						
		int ia = age2530[i2530] - INIT_age; 
		
		for (int ieduccat = 0; ieduccat < 4; ieduccat++){
				
			double n_ieduc = 0.0; 
			int n_ever = 0; 	
			
			for (int iper =0; iper < NUM_smmPerson; iper++){
			
			if (sm_add[ia+1][iper] > 0.01) {
		
				n_ieduc += sm_educcat[i2530][ieduccat][iper]; 
				
				n_ever ++; 
		
			}
			}				
			
			
			double frac_ieduc = 0; 
			if (n_ever >= 1) frac_ieduc = n_ieduc/n_ever; 
			
			if (ieduccat < 2) moments_vec[ivec] = frac_ieduc;
			if (ieduccat < 2) ivec++; 
		
		}
				
	}
	
	
	
			
	int n_everdq_age[30-INIT_age+1]; 
	
	for (int age = INIT_age; age <= 30; age++ ){ 	
		
		int ia = age - INIT_age;
		
		n_everdq_age[ia] = 0; 
		
		for (int iper =0; iper < NUM_smmPerson; iper++){
				
				// include current choice of dq 
				// sm_add0 is an indicator of add=0
			if (sm_add[ia+1][iper] > 0.01){ 
				n_everdq_age[ia] += 1.0;
			} 
		}
		
		// if (age==30) 
		moments_vec[ivec] = 1.0*n_everdq_age[ia]/(1.0*NUM_smmPerson); 
		ivec++; 
	
	}
	
	
	for (int i = numcat_savings-1; i < numcat_savings; i++){
		
		int ia = age_savings[i] - INIT_age; 
		
		// ... savings
		if (age_savings[i] >= 30) moments_vec[ivec] = my_pct(NUM_smmPerson, sm_savings[ia], 50);	
		if (age_savings[i] >= 30) ivec++; 		
					
	}


	const int iage30 = 30 - INIT_age; 

	moments_vec[ivec] = my_gini(NUM_smmPerson, sm_savings[iage30])  ; 	
	ivec++; 


	int iss = 0;  


	int min_ia_logwage = 23-INIT_age; 
	int max_ia_logwage = 34-INIT_age; 
	int n_ols_logwage = 0; 

	iss = 0; 	
	for (int ia = min_ia_logwage; ia <= max_ia_logwage; ia++){
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
			n_ols_logwage += ( sm_earnings[ia][iper] > 1.0e-6 && sm_de[ia][iper] != 1 && sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper]); 
			
			#ifdef DEBUG
				if ( sm_wage[ia][iper] <= 1.0e-4 && sm_work_full[ia][iper] == 1 ) {
					printf("\n\n \t Errors in wage and employment: \t sm_wage[ia][iper] %f,  sm_work_full %d \n exiting.. \n", sm_wage[ia][iper], sm_work_fullpart[ia][iper]); 
					exit(1); 
				}
			#endif 
		}
	}
    
    if (n_ols_logwage < 1) n_ols_logwage = 1;
		
	
	// does not include sd of errors 
	int nx_ols_logwage = NUM_smmcoef_wage - 1; 
	olsY_logwage		= gsl_vector_alloc(n_ols_logwage); 
	olsX_logwage    	= gsl_matrix_alloc(n_ols_logwage, nx_ols_logwage) ;
	olsCoef_logwage	= gsl_vector_alloc(nx_ols_logwage); 	
	gsl_vector_set_zero(olsCoef_logwage); 
	
	iss = 0; 
	for (int ia = min_ia_logwage; ia <= max_ia_logwage; ia++){
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
		if ( sm_earnings[ia][iper] > 1.0e-6 && sm_de[ia][iper] != 1 && sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper]) { 
						
			gsl_vector_set(olsY_logwage, iss, log(sm_earnings[ia][iper])  ); 		
			
			double po_expr = ia+INIT_age - 6.0 - sm_educ[ia][iper]; 
			if (sm_educ[ia][iper]<12) po_expr = ia+INIT_age - 6.0 - 12.0; 
			if (po_expr < 0.0) po_expr = 0.0; 
			
			
			int iss_x = 0; 
			gsl_matrix_set(olsX_logwage, iss, iss_x, po_expr ); 									iss_x++; 		
			gsl_matrix_set(olsX_logwage, iss, iss_x, po_expr*po_expr ); 		iss_x++; 

			gsl_matrix_set(olsX_logwage, iss, iss_x, po_expr * (sm_educ[ia][iper] >= 16) ); 									iss_x++; 		
						
			gsl_matrix_set(olsX_logwage, iss, iss_x, 1.0*(sm_educ[ia][iper] < 12.0 ) * (sm_educ[ia][iper] - 12.0) ); 							iss_x++; 
			gsl_matrix_set(olsX_logwage, iss, iss_x, 1.0*(sm_educ[ia][iper] >= 12.0 ) * (sm_educ[ia][iper] - 12.0) ); 							iss_x++; 
			gsl_matrix_set(olsX_logwage, iss, iss_x, 1.0*(sm_educ[ia][iper] >= 16.0 ) * (sm_educ[ia][iper] - 16.0) ); 							iss_x++; 
			
			gsl_matrix_set(olsX_logwage, iss, iss_x, 1.0*(sm_educ[ia][iper] >= 12 && sm_educ[ia][iper] <= 15) ); 						iss_x++; 
			gsl_matrix_set(olsX_logwage, iss, iss_x, 1.0*(sm_educ[ia][iper] >= 16 ) ); 						iss_x++; 			
			
			gsl_matrix_set(olsX_logwage, iss, iss_x, sm_theta[0][iper]  ); 						iss_x++; 			
			
			gsl_matrix_set(olsX_logwage, iss, iss_x, sm_theta[1][iper]  ); 						iss_x++; 	
			
			gsl_matrix_set(olsX_logwage, iss, iss_x, 1.0 ); 												iss_x++; 	
			
			iss++; 
			
		}
	
	}}

	double chisq_logwage = 0.0; 
	chisq_logwage = olsReg(n_ols_logwage, nx_ols_logwage,   olsY_logwage,    olsX_logwage,   olsCoef_logwage);	
	olsReg(n_ols_logwage, nx_ols_logwage,   olsY_logwage,    olsX_logwage,   olsCoef_logwage);

	
	int min_ia_dedq = 18-INIT_age; 
	int max_ia_dedq = 24-INIT_age; 
	
	int n_ols_dedq = 0;
	for (int ia = min_ia_dedq; ia <= max_ia_dedq; ia++){
    	for (int iper =0; iper < NUM_smmPerson; iper++){

			if ( sm_de[ia-1][iper] == 1 )
			{ n_ols_dedq++; };
	}}
    
    int nx_ols_dedq = 12; 
    gsl_vector * olsY_dedq = gsl_vector_alloc(n_ols_dedq); 
	gsl_matrix * olsX_dedq  = gsl_matrix_alloc(n_ols_dedq, nx_ols_dedq) ;

	gsl_vector * olsCoef_dedq	    = gsl_vector_alloc(nx_ols_dedq); 	    

	iss = 0;  	
	for (int ia = min_ia_dedq; ia <= max_ia_dedq; ia++){
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
			if ( sm_de[ia-1][iper] == 1 )
			{
			
			int iss_x = 0; 
			
			gsl_vector_set(olsY_dedq,    iss, 1.0*sm_dq[ia-1][iper]  ); 		

					
			iss_x = 0; 

			gsl_matrix_set(olsX_dedq,    iss, iss_x, log(sm_cprice[ia-1][iper])  ); 		iss_x++; 

			gsl_matrix_set(olsX_dedq,    iss, iss_x, sm_theta[0][iper]  ); 		iss_x++; 		
			gsl_matrix_set(olsX_dedq,    iss, iss_x, sm_theta[1][iper]  ); 		iss_x++; 	

			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(sm_educ[ia][iper]>=12)  ); 			iss_x++;	
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(sm_educ[ia][iper]>=16)  ); 			iss_x++;	
														
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==18)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==19)  ); 			iss_x++; 			
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==20) ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==21)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==22)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==23)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==24)  ); 			iss_x++; 
			
		
			iss++; 		
			
			}				
	}} 	
	// first stage
	olsReg(n_ols_dedq, nx_ols_dedq,   olsY_dedq,    olsX_dedq,   olsCoef_dedq);
	
	
	iss = 0;  	
	for (int ia = min_ia_dedq; ia <= max_ia_dedq; ia++){
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
			if ( sm_de[ia-1][iper] == 1 ) 
			{
			
			int iss_x = 0; 
			
			gsl_vector_set(olsY_dedq,    iss, 1.0*sm_de[ia][iper]  ); 		
			
			iss_x = 0; 

			double x[] = {log(sm_cprice[ia-1][iper]), sm_theta[0][iper], sm_theta[1][iper], 1.0*(sm_educ[ia][iper]>=12), 1.0*(sm_educ[ia][iper]>=16), 
					1.0*(ia+INIT_age==18), 1.0*(ia+INIT_age==19), 1.0*(ia+INIT_age==20), 1.0*(ia+INIT_age==21), 1.0*(ia+INIT_age==22), 1.0*(ia+INIT_age==23), 1.0*(ia+INIT_age==24)}; 
			
			double dqhat  = 0.0;
			for (int ill=0; ill<nx_ols_dedq; ill++)  dqhat += gsl_vector_get(olsCoef_dedq, ill) * x[ill];  
			
			
			gsl_matrix_set(olsX_dedq,    iss, iss_x, dqhat  ); 							iss_x++; 		// sm_dq_prev[ia][iper]

			gsl_matrix_set(olsX_dedq,    iss, iss_x, sm_theta[0][iper]  ); 		iss_x++; 		
			gsl_matrix_set(olsX_dedq,    iss, iss_x, sm_theta[1][iper]  ); 		iss_x++; 	

			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(sm_educ[ia][iper]>=12)  ); 			iss_x++;	
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(sm_educ[ia][iper]>=16)  ); 			iss_x++;	
														
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==18)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==19)  ); 			iss_x++; 			
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==20) ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==21)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==22)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==23)  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dedq,    iss, iss_x, 1.0*(ia+INIT_age==24)  ); 			iss_x++; 

			iss++; 		

			}				
	}} 	
	
	// second stage
	olsReg(n_ols_dedq, nx_ols_dedq,   olsY_dedq,    olsX_dedq,   olsCoef_dedq);


	
	int min_ia_dqmore = 18-INIT_age; 
	int max_ia_dqmore = 30-INIT_age; 
	
	int n_ols_dqmore = 0;
	for (int ia = min_ia_dqmore; ia <= max_ia_dqmore; ia++){
    	for (int iper =0; iper < NUM_smmPerson; iper++){

			if ( sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper] ){ n_ols_dqmore++; };
	}}
    
    int nx_ols_dqmore_dqlag = 18; 
    gsl_vector * olsY_dqmore_dqlag1 = gsl_vector_alloc(n_ols_dqmore); 
    gsl_vector * olsY_dqmore_dqlag2 = gsl_vector_alloc(n_ols_dqmore); 
	gsl_matrix * olsX_dqmore_dqlag  = gsl_matrix_alloc(n_ols_dqmore, nx_ols_dqmore_dqlag) ;
	gsl_vector * olsCoef_dqmore_dqlag1	    = gsl_vector_alloc(nx_ols_dqmore_dqlag); 	
	gsl_vector * olsCoef_dqmore_dqlag2	    = gsl_vector_alloc(nx_ols_dqmore_dqlag); 	    


	iss = 0;  	
	for (int ia = min_ia_dqmore; ia <= max_ia_dqmore; ia++){
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
			if ( sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper] )
			{
			
			int iss_x = 0; 
			
			gsl_vector_set(olsY_dqmore_dqlag1,    iss, sm_dq[ia-1][iper]  ); 		
			gsl_vector_set(olsY_dqmore_dqlag2,    iss, sm_dq[ia+1][iper]  ); 	
			
			iss_x = 0; 
					
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_educ [ia][iper]  ); 				iss_x++; 

			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_theta[0][iper]  ); 			iss_x++; 		
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_theta[1][iper]  ); 			iss_x++; 	
						
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia][iper])  ); 			iss_x++; 				
			
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia][iper]  ); 			iss_x++; 				

			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia-1][iper]  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia-2][iper]  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia-3][iper]  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia+1][iper]  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia+2][iper]  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, sm_tax_cprice[ia+3][iper]  ); 			iss_x++; 
			
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia-1][iper])  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia-2][iper])  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia-3][iper])  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia+1][iper])  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia+2][iper])  ); 			iss_x++; 				
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, log(sm_cprice[ia+3][iper])  ); 			iss_x++; 
			
			gsl_matrix_set(olsX_dqmore_dqlag,    iss, iss_x, 1.0 ); 									iss_x++;

			iss++; 		

			}				
	}} 	
	
	olsReg(n_ols_dqmore, nx_ols_dqmore_dqlag,   olsY_dqmore_dqlag1,    olsX_dqmore_dqlag,   olsCoef_dqmore_dqlag1);
	olsReg(n_ols_dqmore, nx_ols_dqmore_dqlag,   olsY_dqmore_dqlag2,    olsX_dqmore_dqlag,   olsCoef_dqmore_dqlag2);

	 
	int nx_ols_dqmore = 7; 		
	olsY_dqmore		= gsl_vector_alloc(n_ols_dqmore); 
	olsX_dqmore    	= gsl_matrix_alloc(n_ols_dqmore, nx_ols_dqmore) ;
	olsCoef_dqmore	    = gsl_vector_alloc(nx_ols_dqmore); 	
	gsl_vector_set_zero(olsCoef_dqmore); 
	
	iss = 0;  	
	for (int ia = min_ia_dqmore; ia <= max_ia_dqmore; ia++){
		for (int iper =0; iper < NUM_smmPerson; iper++){
			
			if ( sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper] )
			{
			
			int iss_x = 0; 
			
			gsl_vector_set(olsY_dqmore,    iss, sm_dq[ia][iper]  ); 		
			
			iss_x = 0; 

			double x[] = {1.0*sm_educ [ia][iper] , sm_theta[0][iper], sm_theta[1][iper] , log(sm_cprice[ia][iper]), sm_tax_cprice[ia][iper], 
							sm_tax_cprice[ia-1][iper], sm_tax_cprice[ia-2][iper], sm_tax_cprice[ia-3][iper], sm_tax_cprice[ia+1][iper], sm_tax_cprice[ia+2][iper], sm_tax_cprice[ia+3][iper], 
							log(sm_cprice[ia-1][iper]) , log(sm_cprice[ia-2][iper]), log(sm_cprice[ia-3][iper]), log(sm_cprice[ia+1][iper]), log(sm_cprice[ia+2][iper]), log(sm_cprice[ia+3][iper]), 1.0}; 
			
			double ldqhat  = 0.0;
			for (int ill=0; ill<nx_ols_dqmore_dqlag; ill++)  ldqhat += gsl_vector_get(olsCoef_dqmore_dqlag1, ill) * x[ill];  
			
			
			double l2dqhat = 0.0; 
			for (int ill=0; ill<nx_ols_dqmore_dqlag; ill++)  l2dqhat += gsl_vector_get(olsCoef_dqmore_dqlag2, ill) * x[ill];  		
			
			
			gsl_matrix_set(olsX_dqmore,    iss, iss_x, ldqhat  ); 				iss_x++; 		// sm_dq_prev[ia][iper]
			gsl_matrix_set(olsX_dqmore,    iss, iss_x, l2dqhat ); 					iss_x++; 		// sm_dq_prev[ia-1][iper]
			
								
			gsl_matrix_set(olsX_dqmore,    iss, iss_x, sm_educ [ia][iper]  ); 				iss_x++; 

			gsl_matrix_set(olsX_dqmore,    iss, iss_x, sm_theta[0][iper]  ); 			iss_x++; 		
			gsl_matrix_set(olsX_dqmore,    iss, iss_x, sm_theta[1][iper]  ); 			iss_x++; 	
						
			gsl_matrix_set(olsX_dqmore,    iss, iss_x, log(sm_cprice[ia][iper])  ); 			iss_x++; 				
			

			gsl_matrix_set(olsX_dqmore,    iss, iss_x, 1.0 ); 									iss_x++;

			iss++; 		

			}				
	}} 	
	
	gsl_vector_free( olsY_dqmore_dqlag1 ); 
    gsl_vector_free( olsY_dqmore_dqlag2  ); 
	gsl_matrix_free( olsX_dqmore_dqlag   ) ;
	gsl_vector_free( olsCoef_dqmore_dqlag1	 ); 	
	gsl_vector_free( olsCoef_dqmore_dqlag2	 ); 	    

	olsReg(n_ols_dqmore, nx_ols_dqmore,   olsY_dqmore,    olsX_dqmore,   olsCoef_dqmore);



	
	int n_ols_savings = 0; 
			
	
	for (int iper =0; iper < NUM_smmPerson; iper++){
		for (int i = 1; i < numcat_savings; i++){
		
			int ia = age_savings[i] - INIT_age; 
			
			if ( ia + INIT_age == 30 && sm_savings[ia][iper] >= 100.0 && sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper] ){ n_ols_savings++; };
	}}
    
    if (n_ols_savings<1) n_ols_savings = 1;


	int nx_ols_savings = 3; // NUM_smmcoef_savings; 	
	olsY_savings		= gsl_vector_alloc(n_ols_savings); 
	olsX_savings    	= gsl_matrix_alloc(n_ols_savings, nx_ols_savings) ;
	olsCoef_savings	= gsl_vector_alloc(nx_ols_savings); 
	gsl_vector_set_zero(olsCoef_savings); 
	
	int nx_ols_savings_v2 = 4; 
	olsX_savings_v2    	= gsl_matrix_alloc(n_ols_savings, nx_ols_savings_v2) ;
	olsCoef_savings_v2	= gsl_vector_alloc(nx_ols_savings_v2); 
	gsl_vector_set_zero(olsCoef_savings_v2); 


	
	iss = 0;  	
	for (int iper =0; iper < NUM_smmPerson; iper++){

		for (int i = 1; i < numcat_savings; i++){
		
			int ia = age_savings[i] - INIT_age; 
			
			if ( ia + INIT_age == 30 && sm_savings[ia][iper] >= 100.0 && sm_educ[ia][iper] >= sm_educ[30-INIT_age][iper] ){ 
			
			gsl_vector_set(olsY_savings, iss, log( sm_savings[ia][iper] ) ); 		
			
			int iss_x = 0; 
						
			gsl_matrix_set(olsX_savings, iss, iss_x, sm_theta[0][iper] ); 						iss_x++; 			
			gsl_matrix_set(olsX_savings, iss, iss_x, sm_theta[1][iper] ); 						iss_x++; 
			gsl_matrix_set(olsX_savings, iss, iss_x, 1.0 ); 											iss_x++; 	
			
			
			
			iss_x = 0; 
			gsl_matrix_set(olsX_savings_v2, iss, iss_x, sm_theta[0][iper] ); 						iss_x++; 			
			gsl_matrix_set(olsX_savings_v2, iss, iss_x, sm_theta[1][iper] ); 						iss_x++; 
			gsl_matrix_set(olsX_savings_v2, iss, iss_x, sm_add[ia][iper] ); 											iss_x++; 				
			gsl_matrix_set(olsX_savings_v2, iss, iss_x, 1.0 ); 											iss_x++; 	
			
			
			iss++; 			
			
			}
		
		}
	}	
	
//	olsReg(n_ols_savings, nx_ols_savings,   olsY_savings,    olsX_savings,   olsCoef_savings);
	olsReg(n_ols_savings, nx_ols_savings_v2,   olsY_savings,    olsX_savings_v2,   olsCoef_savings_v2);
	

	
	

		
	for (int iss_x = 0; iss_x < nx_ols_logwage; iss_x++){	
	if (iss_x < nx_ols_logwage-1) {
	moments_vec[ivec] =  gsl_vector_get(olsCoef_logwage, iss_x);
	ivec++;	
	}
	
	}
	moments_vec[ivec] = sqrt(chisq_logwage/n_ols_logwage);
	ivec++;	
	

	for (int iss_x = 0; iss_x < nx_ols_dedq; iss_x++){
	if (iss_x < nx_ols_dedq-9) {
	moments_vec[ivec] = gsl_vector_get(olsCoef_dedq, iss_x);
	ivec++;	
	}
	}
	

	for (int iss_x = 0; iss_x < nx_ols_dqmore; iss_x++){
	if (iss_x < nx_ols_dqmore-1) {
	moments_vec[ivec] = gsl_vector_get(olsCoef_dqmore, iss_x);
	ivec++;	
	}
	}


	
	for (int iss_x = 0; iss_x < nx_ols_savings_v2; iss_x++){
	if (iss_x < nx_ols_savings_v2-1) {
	moments_vec[ivec] = gsl_vector_get(olsCoef_savings_v2, iss_x);
	ivec++;	
	}
	}	

	
	int _SMMnum  = 0;  
	double   _SMMdiff = 0.0;  

	for (int ismm=0; ismm < NUM_ismm; ismm++){

		f[ismm] = moments_vec[ismm] - TARGET_moments[ismm];
        
       
        sumf2     += f[ismm]*f[ismm]*TARGET_wgtmoments[ismm];
		
		        
		SMM_moments[ismm] = moments_vec[ismm]; 
		
		double _absdiff = sqrt(f[ismm]*f[ismm]*TARGET_wgtmoments[ismm]); 
		
		_SMMdiff += _absdiff; 
		_SMMnum  += ( _absdiff >1.0e-6); 
		
	}	
	_SMMdiff /= _SMMnum; 

	#ifdef DEBUG

	printf ("\n \t NUM: %d, ave relative distance %f, \t smm_obj %.5f, %d \n\n", _SMMnum, _SMMdiff, sumf2, NUM_ismm);
		
	if ( ivec != NUM_ismm) {
		printf("\n\n\n !!!!! Error: ivec %d != NUM_ismm %d!!!!!! \n\n exit(1)", ivec, NUM_ismm); 
		exit(1); 
	}
	
		
	#endif 

	
	gsl_vector_free(olsY_logwage ); 
	gsl_matrix_free(olsX_logwage ); 
	gsl_vector_free(olsCoef_logwage ); 
	
	gsl_vector_free( olsY_dedq ); 
	gsl_matrix_free( olsX_dedq   ) ;
	gsl_vector_free( olsCoef_dedq ); 
		
	gsl_vector_free(olsY_savings ); 
	gsl_matrix_free(olsX_savings ); 
	gsl_vector_free(olsCoef_savings ); 
	
//	gsl_vector_free(olsY_savings_v2 ); 
	gsl_matrix_free(olsX_savings_v2 ); 
	gsl_vector_free(olsCoef_savings_v2 ); 
	
	gsl_vector_free(olsY_dqmore ); 
	gsl_matrix_free(olsX_dqmore ); 
	gsl_vector_free(olsCoef_dqmore ); 
		
	
	#ifdef USE_omp
			
	#pragma omp barrier
	#pragma omp flush
			
	#endif	
	
	
	return sumf2; 
	
}


int estimation_nmsimplex(size_t Nx, const gsl_vector *x, gsl_vector *opt_x, gsl_vector *opt_Fval_ss, double init_size, const double size_tolerance, size_t max_iter){
		
	gsl_vector * x_iter = gsl_vector_alloc(Nx);
	
    const gsl_multimin_fminimizer_type *T =  gsl_multimin_fminimizer_nmsimplex;   // gsl_multimin_fminimizer_nmsimplex2;
    gsl_multimin_fminimizer *s = NULL;
    gsl_multimin_function minex_func;
   	
    // Initialize method and iterate
    minex_func.f = &smm_obj_est;
    minex_func.n = Nx;
    minex_func.params = (void *) &Nx; 
	
    // Set initial step sizes according to init_size   
    gsl_vector *ss = gsl_vector_alloc (Nx);
    gsl_vector_set_all (ss, init_size);
	
    s = gsl_multimin_fminimizer_alloc (T, Nx);
    gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

    size_t iter = 0;    int status = 0;    double size = init_size;		
    do {
	
		#ifdef USE_omp
		#pragma omp barrier
		#pragma omp flush
		#endif 
		
	 	iter++; 
	 
	 	status = gsl_multimin_fminimizer_iterate(s);
	
		if (status) break;
   	
        size   = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, size_tolerance);
				         
		if (status == GSL_SUCCESS) printf ("converged to minimum at\n");
		
		// if ( iter%1==0 )
		{ 
			printf ("%zu \t f = %f, \t size = %f \n", iter, s->fval, size);
			
			for (size_t i_s =0; i_s<Nx; i_s++){
				printf(" \t x [%3zu] = %5.6f ; \n", i_s, gsl_vector_get (s->x, i_s));
				gsl_vector_set(x_iter, i_s, gsl_vector_get (s->x, i_s)); 
			}
			
			getpar(x_iter); printpar("estpar_iter.txt"); 
			printvec("moments_vec_iter.txt", SMM_moments, NUM_ismm);	
						
			if ( mpi_id == 0) { getpar(x_iter); printpar("estpar_iter_check.txt"); }
			
			printvec("estvec_iter.txt", x_iter, Nvar);
			
		}
		
	} while (status == GSL_CONTINUE && iter < max_iter);

    
	// store the optimal points;
	for (size_t i=0; i<Nx; i++) gsl_vector_set(opt_x, i, gsl_vector_get (s->x, i)); 	
	
	// store the function value at optimal point and the step size
	gsl_vector_set(opt_Fval_ss, 0, (s->fval)); 
	gsl_vector_set(opt_Fval_ss, 1, gsl_multimin_fminimizer_size (s) ); 
	
	printf("\n Write the estimated parameter values to file estpar_iter.txt \n");
	getpar(opt_x); 
	printpar("estpar_iter.txt");
	
    gsl_vector_free(ss);  gsl_multimin_fminimizer_free (s); 
	gsl_vector_free(x_iter); 
	
    return status;

}


void estiamte(gsl_vector *vecx){

	// !!!!!!!!!!! PARAMETER ESTIMATION STARTS !!!!!!!!!!	
			
	// Start estimation: 
	printf ("\n\n Estimation Starts!!! \n");
	
	time_t start0,end0;	double timediff; 
	
	time (&start0);
	
	double temp_dis = smm_obj(vecx,vecx); 
	
	time (&end0);	timediff = difftime (end0,start0);
	
	printf ("\n\n It took you %f minutes to evaluate moment conditions, Initial distance = %f \n", timediff/60, temp_dis);
	
	
	time_t start_t, end_t;	double diff_t; 
	time (&start_t);
	gsl_vector *opt_vecx 	= gsl_vector_alloc(Nvar);	
	gsl_vector *opt_f 	= gsl_vector_alloc(2);
	estimation_nmsimplex(Nvar, vecx, opt_vecx, opt_f, Est_stepsize, Est_tol, Est_maxiter);
	time (&end_t);	diff_t = difftime (end_t, start_t);
	printf ("\n It took you %f hours to finish Estimation: f %f, size %f \n", diff_t/3600, gsl_vector_get(opt_f,0), gsl_vector_get(opt_f, 1) );
	
	// !!!!!!!!!!! PARAMETER ESTIMATION COMPLETES !!!!!!!!!!	
	
	double opt_dis = smm_obj(opt_vecx, opt_vecx); 
	printf("\n check opt_dis: %f", opt_dis); 
	
    gsl_vector_memcpy(vecx, opt_vecx);
    
	gsl_vector_free(opt_vecx); gsl_vector_free(opt_f); 
		
}
